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ABSTRACT 

The quantitative morphological classification of distant galaxies is essential to the understanding of the 
evolution of galaxies over the history of the Universe. This paper presents Hubble Space Telescope WFPC2 
F606W and F81AW photometric structural parameters for 7450 galaxies in the "Groth Strip." These 
parameters are based on a two-dimensional bulge+disk surface brightness model and were obtained using 
an automated reduction and analysis pipeline described in detail here. A first set of fits was performed 
separately in each bandpass, and a second set of fits was performed simultaneously on both bandpasses. 
The information produced by these two types of fits can be used to explore different science goals. 
Systematic and random fitting errors in all structural parameters as well as bulge and disk colors are 
carefully characterized through extensive sets of simulations. The results of these simulations are given 
in catalogs similar to the real science catalogs so that both real and simulated measurements can be 
sampled according to the same selection criteria to show biases and errors in the science data subset of 
interest. The effects of asymmetric structures on the recovered bulge+disk fitting parameters are also 
explored through simulations. The full multidimensional photometric survey selection function of the 
Groth Strip is also computed. This selection function, coupled to bias maps from simulations, provides 
a complete and objective reproduction of the observational limits, and these limits can be applied to 
theoretical predictions from galaxy evolution models for direct comparisons with the data. 

Subject headings: galaxies : fundamental parameters, galaxies : evolution 



1. INTRODUCTION 

The visual classification of galaxies has a venerable tradi- 
tion in optical astronomy starting with the introduction of 
Hubble's famous "Tuning Fork" diagram (Hubble 1936). 
Despite the fact that others have extended Hubble's origi- 
nal classification mainly to deal with the diversity of struc- 
tures in later-type galaxies (de Vaucouleurs 1959; van den 
Bergh 1960a, b, 1976; Morgan 1970), the nomenclature of 
Hubble still very much pervades the language of today's 
galaxy morphology literature. This longevity is a trib- 
ute to Hubble's seminal work. However, visual classifica- 
tion has serious weaknesses. First and foremost, it is a 
subjective process. Although visual classification experts 
can agree to within two revised Hubble types (Nairn et al. 
1995), it is very difficult for non-experts to produce reliable 
visual classifications without years of hard- won experience. 
Second, it is also unclear how useful visual classification is 
with respect to high-redshift galaxies. Limited spatial res- 
olution means that larger and larger internal galaxy struc- 



tures such as spiral arms and tidal tails get progressively 
smoothed out with increasing redshift, and this smoothing 
can introduce significant classification biases. 

No visual or quantitative classification system is perfect. 
However, provided a given quantitative system is clearly 
defined by its proponents, it is readily reproducible in its 
successes and failures by others. This is the fundamental 
advantage of the quantitative approach. Moreover, sys- 
tematic and random errors of quantitative classifiers can 
be carefully characterized through extensive sets of galaxy 
image simulations covering a wide range of structural pa- 
rameters. This is another important advantage. It should 
be emphasized that quantitative classification is meaning- 
less without three important elements: the measurements 
themselves, the simulations and the galaxy selection func- 
tion. This selection function is critical to relate observed 
structural parameter distributions to predictions from the- 
oretical models. 

A number of quantitative classifiers have been devel- 
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oped/extended over the years to probe the structure of 
high-redshift galaxies. These classifiers can be either para- 
metric (model-based) or non-parametric. Non-parametric 
classifiers include the C — A system (Watanabe et al. 1985; 
Abraham et al. 1994, 1996; Wu 1999; Conselice, Bershady, 
& Jangrcn 2000) , artificial neural nets trained from visual 
classification sets (Odewahn 1995; Odcwahn et al. 1996), 
and self-organizing maps (Nairn et al. 1997). Paramet- 
ric classifiers include radial multi-gaussian deconvolution 
(Bendinelli 1991; Fasano et al. 1998) and bulge+disk de- 
composition (Schade et al. 1995, 1996; Ratnatunga et al. 

1999) . The latter is popular for three reasons: (1) It is 
rooted in the very first studies of the functional form of 
galaxy radial surface brightness profiles (de Vaucouleurs 
1948, 1959); (2) It provides a "comfortable" mental pic- 
ture of the overall structure of a distant galaxy, i. e., it 
is conceptually simple to relate a quantitative measure- 
ment of galaxy type such as bulge-to-total light ratio to 
the familiar Hubble types; and (3) Photometric entities 
such as bulges and disks have distinct dynamical counter- 
parts although this correspondence does not always hold 
(see Section 5.2 for more details). 

This work uses GIM2D (Galaxy IMage 2D), a 2D de- 
composition fitting program (Simard 1998, and this pa- 
per), to measure the structural parameters of galaxies 
in the "Groth Strip" (Groth et al. 1994; Rhodes et al. 

2000) . GIM2D is an IRAF 9 /SPP package written to per- 
form detailed bulge+disk surface brightness profile decom- 
positions of low signal-to- noise (S/N) images of distant 
galaxies in a fully automated way. GIM2D takes an HST 
or ground-based science image and its source catalog, per- 
forms 2D profile fits on each source and produces model- 
subtracted images as well as a full catalog of structural pa- 
rameters. The bulge+disk model adopted here is not fun- 
damentally new, but GIM2D offers an independent check 
of other galaxy classification works by including a set of 
extended features (Sersic bulge profile, a comprehensive 
but by all means not exhaustive set of image asymmetry 
indices, three different fitting methods) and a different fit- 
ting algorithm. GIM2D has already been used in a variety 
of HST and ground-based distant galaxy studies: the op- 
tical structure of intermediate redshift compact narrow- 
emission line galaxies (Guzman et al. 1998), the quan- 
titative morphology of Hubble Deep Field North galax- 
ies (Marleau & Simard 1998), the NICMOS structure of 
a spiral galaxy lens at z = 0.4 (Mailer et al. 2000), the 
luminosity-size relation of field disk galaxies from z = 0.1 
to z = 1.1 (Simard et al. 1999), the number density and 
luminosity function of E/SO galaxies to z < 1 (Im et 
al. 2001), the Fundamental Plane of field absorption- line 
galaxies out to z <~ 1 (Gcbhardt et al. 2002), tests of hi- 
erarchical galaxy evolution models (Simard et al. 2002), 
the colors of luminous bulges at z ~ 1 (Koo et al. 2002), 
the galaxy populations of poor, X-ray selected groups of 
galaxies (Tran et al. 2001) and of high and low X-ray lu- 
minosity galaxy clusters (Balogh et al. 2002), and the size 
evolution of high-redshift brightest cluster galaxies (Nel- 
son et al. 2002). 

This paper describes in detail the structural parameter 
analysis of galaxies in the Groth Strip from calibrated HST 

9 IRAF is distributed by the National Optical Astronomy Observato 
Astronomy, Inc., under cooperative agreement with the National Sci 



archival images to final, parameter-rich structural catalogs 
for the entire Strip. It is a companion paper to Vogt et al. 
(2002) and Phillips et al. (2002) in which the spectroscopic 
Keck/Low Resolution Imaging Spectrograph survey of the 
Groth Strip undertaken by the Deep Extragalactic Evo- 
lutionary Probe (DEEP) team is described. Cosmological 
parameters adopted throughout this paper are H a = 70 
km s- 1 Mpc" 1 , and (Ct m , Ct A , Ct k ) = (0.3, 0.7, 0.0). 

2. OVERVIEW 

Quantitative galaxy morphology classifiers must include 
three vital elements in order to be scientifically useful: 
the structural measurements themselves, extensive simu- 
lations, and a detailed survey selection function. 

Every structural catalog should have a companion cata- 
log of simulations from which the simulated counterparts 
of real data plots can be extracted. The ability to clearly 
visualize the systematic and random errors of a set of ob- 
served structural parameters is very helpful to quickly as- 
sess what science goals can be effectively pursued with 
the measurements. Morever, the systematic and random 
errors as characterized through these simulations can be 
applied to the theoretical predictions from galaxy evolu- 
tion models to "scatter" the models in the same way as 
the observational errors would scatter signals from the real 
Universe. 

The survey selection function is constructed by insert- 
ing galaxy images with a wide range of input structural 
parameters into real images and calculating the success 
rate of the detection algorithm as a function of those pa- 
rameters. This selection function serves two purposes: (1) 
It can be used to insure that studies of distant galaxies 
spanning a large range of redshifts treat galaxy samples 
at different redshifts in the same way (e.g., Simard et al. 
1999). (2) The selection function can be used to "observe" 
theoretical models for a direct comparison with the data 
(e.g., Simard et al. 2002). Together with the simulations, 
the selection function provides the best possible way to 
reproduce the biases of the observational strategy. 

Current and planned future on-line archival data sys- 
tems contain (or will contain) prodigious amounts of data. 
Mining these systems in as automated a mode as possible 
holds the promise of fantastic scientific returns on prob- 
lems requiring very large, statistically well-defined samples 
such as the evolution of galaxies. The automated GIM2D 
pipeline described in the next sections was designed to pro- 
duce the above three elements with data mining of large 
datasets in mind. 

The GIM2D pipeline includes the following steps: (1) 
Pre-processing, cosmic ray (CR) rejection and data quality 
mask creation (Section 3), (2) Source detection, deblcnd- 
ing and extraction (Section 4), (3) Bulge+disk decompo- 
sitions of galaxy images (Section 5) with a Point-Spread- 
Function (PSF) for each object (Section 5.1), (4) Creation 
of residual images and computation of residual image in- 
dices (Section 5.6), (5) Construction of measured struc- 
tural parameter catalogs (Section 6), (6) Construction of 
associated catalogs of simulations for mapping systematic 
biases and random errors (Section 7), and (7) Construction 
of the survey selection function (Section 8). 

., which arc operated by the Association of Universities for Research in 
c Foundation. 
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3. PRE-PROCESSING OF HST IMAGES 

The GIM2D pipeline starts with the retrieval and pre- 
processing of HST archival images. Archival images are 
calibrated "on-the-fly" at retrieval time. The retrieved 
images are then sent through a set of IRAF scripts that 
removes cosmic rays and combine all the images taken in 
the same bandpass and at the same position on the sky to- 
gether. The pipeline does not mosaic the three WFC chips 
together. Mosaicing (e. g., IRAF/STSDAS/WMOSAIC) 
destroys the uncorrelated noise characteristics of the orig- 
inal images by applying flux interpolations and geometric 
corrections. The noise characteristics of the science im- 
ages must be preserved as much as possible through the 
whole pipeline so that it is possible to provide a realis- 
tic noise estimate to the likelihood function during the 
galaxy image fitting process. Each WFPC2 pointing is 
therefore split into three different images, one for each 
Wide-Field chip, and each chip is then processed inde- 
pendently through the rest of the pipeline. Each WFC 
chip covers a 80" x 80" field of view with a pixel scale of 
O'T/pixel. The Planetary Camera images were not ana- 
lyzed since all the DEEP/GSS Keck spectroscopic targets 
were selected solely from the WFC images for the sake of 
sample homogeneity. 

3.1. HST/WFPC2 Image Datasets 

The HST images are from two surveys, collectively 
dubbed the "Groth Survey Strip" (GSS), taken under 
HST programs GTO 5090 (PI: Groth) and GTO 5109 (PI: 
Westphal). The GSS consists of 28 overlapping subfields 
taken with the HST Wide Field and Planetary Camera 2 
(WFPC2) and forms a "chevron strip" oriented NE to SW 
at roughly 1417+52 at Galactic latitude b ~ 60°. Each of 
27 subfields has exposures of 2800 seconds (4x 700 sec- 
onds) in the broad V filter (F6061V) and 4400 seconds 
(4 x 1100 seconds) in the broad I filter (F814W). The 
28th field is the Westphal Deep Survey Field 2 (J2000 
1417.5+52.5), with total exposures of 24,400 seconds in 
V and 25,200 seconds in I. The images were recalibrated 
"on-the-fly" through the Canadian Astronomy Data Cen- 
ter (CADC) standard pipeline (Pirenne et al. 1998). The 
datasets analyzed in this paper are listed in Table 1. Field 
7 is the Westphal Deep Survey Field. 

3.2. WFPC2 Data Quality Frames 

Archival data produced through on-the-fly recalibration 
come in two parts: the science frames and the data qual- 
ity frames (DQF). Each science frame comes with a data 
quality frame. This frame flags pixels in the science frame 
that have been corrupted for one reason or another (see 
Section 26.2.2 of HST Data Handbook Version 3.0 for DQF 
flag values). A DQF flag value of zero is used by the HST 
pipeline to flag good pixels. The GIM2D pipeline goes 
through the data quality frames of all science images in 
a given image stack 10 and identifies all the pixels with 
no good pixel values in the stack. The locations of these 
pixels are recorded and applied to the SExtractor segmen- 
tation image (Section 4.2) to make sure bad pixels are not 
included in the surface brightness fits. 

3.3. Cosmic Ray Rejection and Image Combination 

10 An image stack is denned here as a set of consecutive exposures taken 



The IRAF/STSDAS task CRREJ was used to reject cos- 
mic rays from GSS WFPC2 image stacks. It combines a 
stack of consecutive exposures while eliminating cosmic 
rays and scaling the remaining pixels to the total expo- 
sure time of the stack. Cosmic rays are rejected through a 
series of sigma-clipping iterations. The number and sigma 
thresholds of these iterations are specified using the pa- 
rameter SIGMAS. SIGMAS must be chosen carefully If 
SIGMAS is too high, too many cosmic rays will be missed. 
On the other hand, if SIGMAS is too low, the CR rejec- 
tion algorithm will start eating away at the noise in the 
background pixels. Unfortunately, there is no prescribed 
way of determining SIGMAS. The approach adopted here 
involved first the creation of a combined image with such 
high SIGMAS values that no pixels were rejected by CR- 
REJ. The resulting "HIGH" image showed all the cosmic 
rays that hit all the science images in the image stack. 
However, there were also plenty of untouched background 
pixels that showed what the background in the combined 
image should look like. CRREJ-combincd images were 
then created with progressively lower values of SIGMAS, 
and blinked against the "HIGH" image. SIGMAS values 
were lowered until the noise in the untouched background 
areas of these images clearly showed they were being mod- 
ified by the sigma-clipping with respect to the same areas 
in the' "HIGH" image. The best value of SIGMAS was 
found to be "6,4" i. e., CRREJ was instructed to perform 
two sigma-clipping iterations, the first one at the 6cr-level 
and the second one at the 4<r-level. This value was then 
automatically applied in the combinations of all the GSS 
stacks. Even with the best SIGMAS cuts, a number of low- 
energy cosmic rays will be left in the final combined image. 
If these low-energy cosmic rays significantly changed the 
background noise properties of the final combined images, 
the background pixel histograms of these images should 
show deviations from a Gaussian distribution in the form 
of high flux tails. No such deviations were ever observed in 
the histograms inspected at different background locations 
in the final images. 

3.4. WFPC2/GSS Photometric Zero-points 

All total F814W fluxes F 814 (galaxy, bulge or disk) will 
be converted in this paper to F814TV magnitudes on the 
Vega system (denoted Isi4 or simply I hereafter) using the 
equation: 

1814 = -2.5 log 10 (F 8 i 4 /t 8 i 4 ) + C 8 14, (1) 

where C 8 i 4 = 21.65 (May 1997 WFPC2 SYNPHOT up- 
date). 

Similarly, all total F606W fluxes F@o6 (galaxy, bulge or 
disk) were converted to F606IV magnitudes on the Vega 
system (denoted Vqog or simply V hereafter) using the 
equation: 

V 60 6 = -2.51og 10 (f 6 o 6 /t606) + Ceoe, (2) 

where C 606 = 22.91 (May 1997 WFPC2 SYNPHOT up- 
date). The total exposure time i 8 i 4 was 4400 seconds in 
the F814FV filter, and the total exposure time teo6 was 

at the same location on the sky and through the same filter. 
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2800 seconds in the F606IU filter in all GSS fields except 
Field 7 for which t S u = 25200 seconds and t 60e = 24400 
seconds. 

4. SOURCE DETECTION, DEBLENDING AND EXTRACTION 

To proceed with the fitting of galaxy images, GIM2D 
needs a catalog of sources for each image to be analyzed. 
For each source the catalog must include a x-y centroid po- 
sition, an initial estimate of the local sky background level 
and the isophotal area of the object in pixels above the 
detection threshold. GIM2D also needs a segmentation or 
mask image in which pixels belonging to the same object 
are all assigned the same flag value and sky background 
pixels are flagged by zeroes. The source catalogs and seg- 
mentation images for the Groth Strip were created using 
the SExtractor ("Source Extractor") galaxy photometry 
package version 1.0a (Bertin & Arnouts 1996). 

4.1. Detection Parameters 

The SExtractor source detection was run on the 
CRREJ-combined /§14 GSS images. The detection thresh- 
old was 1.5-abkg, and the required minimum object area 
above that threshold was 10 pixels. The convolution ker- 
nel was a 3x3 Gaussian kernel with a FWHM of 1.5 pixels. 
These detection parameters do not have to be particularly 
fine-tuned to extract the faintest possible sources from the 
GSS images since the faintest magnitude at which reli- 
able bulge+disk decompositions can be performed is well 
above the magnitude limits of the SExtractor source cat- 
alog. No star/galaxy separation was attempted. Every 
source was fitted with GIM2D. Unresolved sources such as 
stars could easily be identified as GIM2D output models 
with zero half-light radius. 

4.2. Deblending 

As SExtractor performs source detection and photom- 
etry, it is able to dcblcnd sources using flux multi- 
thresholding. This deblending technique works well in the 
presence of saddle points in the light profiles between ob- 
jects. Each SExtractor pre-deblending "object" consists 
of all the pixels above the detection threshold that are 
spatially connected to one another. This group of pixels 
may or may not include several real objects. For each 
"object," the multi-thresholding algorithm goes through 
its connected pixels and rethrcsholds them at N levels 
(N = 32 in the current analysis) between the l.5o~bk g 
isophote and the peak pixel value in the "object" to build 
a two-dimensional flux tree of the "object." The algorithm 
then goes down the threshold levels, and it looks at each 
branch in the tree to see if the flux contained in that branch 
above the threshold level is a fraction f sep or greater of 
the total flux in the "object." If a given branch meets this 
condition, it is then treated as a separate object, and the 
separation threshold for that object is set to the thresh- 
old at which the split occurred. The multi-thresholding 
algorithm assigns the pixels between two adjacent objects 
and below the separation threshold based on a probability 
calculated from bivariate Gaussian fits to the two objects. 
No assumption is made regarding the shape of the objects 
in this statistical deblending technique. 

The fraction f sep is set by the SExtractor input param- 
eter DEBLEND_MINCONT. A value of 0.00075 was used 



for the SExtractor GSS source catalogs. This value is sub- 
jective, and it was found through visual inspection of sev- 
eral GSS fields to provide good object separation. Even 
though the value of DEBLEND_MINCONT was deter- 
mined subjectively, it provides an unequivocal definition 
of an object in the GSS catalogs presented in this pa- 
per. It was only determined once, and the same value 
of DEBLEND.MINCONT was consistently used for all 
GSS fields as well as for all GSS simulations. 

Figure 1 shows a typical section of the Groth Strip taken 
through the F8UW filter (DEEP/GSS Field ID 8/WFC 
Chip 3), and Figure 2 shows the corresponding SExtractor 
segmentation image. 

4.3. Thumbnail Image Extraction 

GIM2D disk+bulge decompositions are performed on 
thumbnail (or "postage stamp" ) images extracted around 
the objects detected by SExtractor rather than on the en- 
tire science image itself. Thumbnail images are preferable 
for two reasons: (1) Thumbnail images reduce the mem- 
ory and I/O footprints of the program so it can be used 
in background mode on many computers without signifi- 
cantly impacting their other users. (2) Many CPUs can 
work on the same science image at the same time as they 
work on different thumbnail images. This is extremely use- 
ful since the same list of thumbnail images can be sent to 
all available CPUs on a network/cluster, and all CPUs will 
work down the same master list without interfering with 
one another. There is no limit on the number of CPUs 
that can be simultaneously harnessed for a given master 
list. 

GIM2D will extract two or three thumbnail images for 
each object in the SExtractor catalog. The area of an ob- 
ject's thumbnail images is given by the isophotal area of 
the object. Here, all thumbnails were chosen to have an 
area 20 times larger than the l.b-o~bkg isophotal area. The 
first thumbnail is extracted from the science image itself, 
and the local background calculated by SExtractor is sub- 
tracted from it so that it should have a background mean 
level close to zero. The second thumbnail is extracted from 
the SExtractor segmentation image. This segmentation 
thumbnail is modified so that bad pixels identified by the 
DQF frames (see section 3.2) will be excluded from the 
surface brightness fits. Some HST image datasets, most 
notably NICMOS datasets, include a "sigma" image giv- 
ing the expected background + Poisson noise at each pixel. 
If such a sigma image is available, GIM2D will automati- 
cally extract the third thumbnail image from it. 

5. SURFACE BRIGHTNESS FITS 

5.1. Point-Spread-Functions 

GIM2D accepts four kinds of Point-Spread-Functions 
(PSFs): a 2D gaussian PSF, a delta function PSF, a user- 
given PSF or a TinyTim PSF. GIM2D automatically nor- 
malizes the total flux in all input PSFs to 1.0 to ensure 
that this step has been performed. The 2D gaussian PSF 
is generated by GIM2D with the seeing FWHM specified 
in the GIM2D parameter file. If the delta function PSF 
option is selected, no PSF-convolution is performed on the 
galaxy image models. PSF effects should not be important 
for structures spanning a large number of resolution ele- 
ments. Fits without PSF convolution require considerably 
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less computation time, a definite advantage for very large 
objects. The third kind of PSF is usually an image given 
by the user. For example, this image could be a bright star 
extracted directly from the image or could be created from 
a set of PSF stars using a stellar photometry program such 
as DAOPHOT (Stetson 1987). This option is particularly 
useful for ground-based studies in which sufficiently sam- 
pled PSF stars are easily found all over the science frames. 
TinyTim PSFs are generated by the Space Telescope Sci- 
ence Institute package TINY TIM (Krist 1993) version 
4.4. TinyTim PSFs are available for all major imaging in- 
struments onboard HST including WPFC2 and NICMOS. 
WFPC2 TinyTim PSFs were used for all (small and large) 
objects in the GSS analysis here. 

The shape of the HST/WFPC2 PSF varies significantly 
as a function of position. Therefore, F814W and F606W 
PSFs were generated with TinyTim for each GSS object 
analyzed here. All GSS PSFs were oversampled by a factor 
of 5 and were 2 '.'4 on a side. No telescope jitter was added. 
The pointing stability of HST in fine lock mode is typi- 
cally better than 5 mas RMS, and any value under 7 mas 
is not noticeable (Krist 1993). WFPC2 charge- coupled 
devices (as any other CCDs) suffer from charge diffusion. 
Charge diffusion contributes a small amount of blurring 
to the images, and TinyTim does include a WFPC2 sub- 
pixel response function to take charge diffusion into ac- 
count. However, it is very important to mention that Tiny- 
Tim does not automatically apply this subpixcl response 
function to oversampled PSFs. GIM2D automatically con- 
volves the oversampled PSF-convolved galaxy model with 
the subpixel response function kernel given in Krist (1993) 
after the model has been rebinned to the WFC detector 
resolution. 

5.2. Galaxy Image Model 

The bulge+disk model used in GIM2D and other works 
is obviously a simple approximation. After all, many real 
galaxies will exhibit more than two structural components 
such as nuclear sources, bars, spiral arms and HII regions. 
Even in the presence of only a bulge and a disk, the el- 
lipticity and/or the position angles of these components 
might be functions of galactocentric distance. However, 
each new structural component or new layer of complex- 
ity added to the model comes with additional parameters 
that stretch the amount of information that can be rea- 
sonably extracted from small, low signal-to-noise images 
of distant galaxies. The simple bulge+disk model is a 
trade-off between a reasonable number of fitting param- 
eters and a meaningful decomposition of galaxy images. 
Despite its relative simplicity, careful analysis of the pa- 
rameters of the bulge+disk model can yield useful infor- 
mation regarding those higher layers of complexity. For 
example, a very large deviation between the position an- 
gles of the bulge and of the disk is usually a strong in- 
dication of the presence of a bar, so barred galaxies can 
be easily identified without individually looking at all the 
galaxies in the sample. Histograms of A</> = \<pb — <pd\ 
are indeed a powerful way of finding barred galaxies. The 
shape of such histograms typically shows a well defined 
peak at A</> = with a tail of outliers that turn out to be 
the barred galaxies. The bulge+disk model is also a good 
way to study the morphology of quasar host galaxies. In 



the presence of a strong, unresolved, central source, the 
radius of one of the model components will shrink to zero 
to accommodate the source, and the other component will 
shape itself after the host galaxy. So, if quasar nuclei are 
found predominantly in elliptical galaxies, the radius of 
the disk component would be the one to converge to zero. 
Differences in the F814W and FQ06W centroid can also 
provide information on the irregularity of galaxy shape. 
All these examples illustrate the fact that the apparent 
simplicity of the bulge+disk model should not mask the 
richness of information provided by the inter-comparison 
of its full set of parameters. 

Nonetheless, the bulge+disk model used here (and in 
other works) also has important limitations: it has a 
unique center for the whole galaxy, and the intrinsic (i. e., 
before PSF convolution) ellipticity and position angle of 
each component do not change as a function of radius. 
This is in contrast to techniques such as isophotal ellipse 
fitting in which it is customary to let ellipticity, position 
angle and centroid vary from one ellipse to the next. The 
limitations of the bulge+disk model can yield results that 
are unexpected (or unsuspected!) at first glance. Consider 
a purely elliptical galaxy with no disk component whatso- 
ever but with varying ellipticity and position angle with 
radius. A bulge+disk fit to such an object may converge 
to a model with a significant disk component (B/T < 1). 
The additional degrees of freedom provided by the disk 
component are used by the model to compensate for the 
varying ellipticity and position angle. The position angle 
difference between the bulge and disk components will thus 
reflect the position angle difference between the inner and 
outer isophotes of the galaxy, and the inclination angle of 
the disk component will reflect the ellipticity of the outer 
isophotes. If mergers in high-density environments such 
as galaxy clusters induce such changes in ellipticity and 
position angle with galactocentric radius, then the same 
range of bulge fraction may not select the same galaxies 
in clusters and in the field. In addition to possible de- 
viations from a pure deVaucouleurs law, this behavior of 
the bulge+disk model explains why, for example, brightest 
cluster galaxies can often have bulge fractions as low as 0.4 
(Nelson et al. 2002). Finally, a single center for the model 
makes it impossible to detect differences in bulge and disk 
centroids should they manifest themselves in some stages 
of the evolution of galaxies. Unfortunately, the above lim- 
itations cannot be adequately dealt with by adding more 
fitting parameters since the basic bulgc+ disk model al- 
ready has enough parameters to make the search for the 
best-fit solution arduous. 

It should be kept in mind that the conventional 
"bulge/disk" nomenclature does not say anything about 
the internal kinematics of the components. The pres- 
ence of a "disk" component does not necessarily imply 
the presence of an actual disk since many dynamically hot 
systems also have simple exponential profiles of the form 
given by Equation 5 below (Lin & Faber 1983; Kormendy 
1985). Likewise a "bulge" may represent a brightened cen- 
ter due to a starburst rather than a genuine dynamically 
hot spheroid. To avoid any confusion between photometric 
structures and internal dynamics entities, the names "pho- 
tobulge" (for "photometric bulge") and "photodisk" (for 
"photometric disk" ) will be used hereafter to refer to the 
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r 1 /™ and exponential components of galaxy light profiles 
respectively. 

The 2D galaxy model used by G1M2D has a maximum of 
twelve parameters: the total flux F in data units (DU), the 
bulge fraction B/T (= for pure photodisk systems), the 
photobulge semi-major axis effective radius r e , the pho- 
tobulge ellipticity e (e = 1 — b/a, b = semi-minor axis, 
a = semi-major axis), the photobulge position angle of 
the major axis (fit, on the image (clockwise, y-axis = 0), 
the photodisk semi-major axis exponential scale length 
(also denoted h in the literature), the photodisk inclina- 
tion i (face-on = 0), the photodisk position angle <frd on 
the image, the subpixel dx and dy offsets of the model 
center with respect to the thumbnail science image cen- 
ter, the background residual level db, and the Sersic index 
n. Twelve parameters is a maximum since one or more 
parameters can be frozen to their initial values if neces- 
sary depending on the scientific goals being pursued. The 
position angles (f>b and <f>d were not forced to be equal for 
two reasons: (1) a large difference between these position 
angles is a signature of barred galaxies, and (2) some ob- 
served galaxies do have bona fide photobulges that are not 
quite aligned with the photodisk position angle. 

Two types of radii are in use in the literature: geomet- 
ric mean (also known as "circular" or "equivalent") radii 
and semi-major axis radii. The same name i. e., "effec- 
tive" is used for both types of radii, and confusion arises 
when the type of radius used for a given sample is not 
clearly specified. As noted above, the photobulge effec- 
tive radius in the GIM2D image model is calculated with 
respect to the semi-major axis of the photobulge compo- 
nent. Since the geometric mean (or "circular") photob- 
ulge effective radius r e ^ circ is given by Vafr, it is related to 
the semi-major axis radius r e . sma by the simple relation 
fe.drc = f e , sma \/l — e. Similarly, the GIM2D semi-major 
axis photodisk scale length rd lS ma is related to the cir- 
cular photodisk scale length by r ejCirc = r e ^ sma ^l - e d isk 
where edisk is the ellipticity of the photodisk. Given the 
thinness of galaxy disks and the great distances of high- 
rcdshift galaxies, the GIM2D photodisk model does not 
need to include a vertical scale height, The inclination an- 
gle of this infinitely thin photodisk is calculated from the 
measured photodisk ellipticity as i — arccos(v / l — edisk) 
assuming that face-on photodisks are perfectly circular. 
This relation between photodisk ellipticity and inclination 
is different from the prescriptions used locally in Tully- 
Fisher work (e.g., Courteau 1997) which must account for 
the effect of disk scale height on the observed ellipticity of 
highly-inclined spiral galaxies. 

The first component ("photobulge") of the 2D surface 
brightness model used by GIM2D to fit galaxy images is a 
Sersic profile of the form: 

E(r) = ^expi-kKr/r^ 1 - 1]}, (3) 
where £(r) is the surface brightness at radius r (Sersic 
1968). The parameter k is set equal to 1.9992n-0.3271 
so that r e remains the projected radius enclosing half of 
the light in this component (Capaccioli 1989). The clas- 
sical de Vaucouleurs profile has the special value n = 4, 
and this value was chosen for the current analysis. This 
choice was motivated by studies of bulge profiles in lo- 
cal galaxies. Locally, there is evidence that the bulges 
of late-type spiral galaxies may be better fitted by an n 



= 1 profile, whereas bright ellipticals and the bulges of 
early-type spiral galaxies follow an n — 4 profile (de Jong 
1996; Courteau et al. 1996; Andredakis 1998). Local late- 
type galaxies with n = 1 bulges have B/T < 0.1 (de Jong 
1996). Since such bulges contain only 10% of the total 
galaxy light, low signal-to-noise measurements of late- type 
high-redshift galaxies make it very difficult to determine 
the Sersic index. On the other hand, n is more important 
for bulge-dominated galaxies, and n — 4 is the expected 
value based on local early- type galaxies. Knowing that 
bright ellipticals and the bulges of early-type spirals are 
well-fitted by a de Vaucouleurs profile, an = 4 bulge pro- 
file was therefore adopted as the canonical bulge fitting 
model here for the sake of continuity across the full range 
of morphological types. The total flux in the Sersic photo- 
bulge component is calculated by integrating Equation 3 
from r = to infinity to obtain: 

^photobulge = 27me fc fc- 2n r e 2 r(2n)£ e (4) 

where F is the gamma function. For n = 4, -Fphotobuige = 
7.2147rrg£ e . The bulge component was collapsed to a 
point (zero radius) source anytime its effective radius was 
less than 0:'02 (0.2 WFC pixel) when the bulge+disk 
model fits were performed. 

The second component ("photodisk") of the GIM2D 
model is a simple exponential profile of the form: 

E(r) = Y. exp(-r/r d ). (5) 

So is the face-on central surface brightness. The photodisk 
is assumed to be infinitely thin. The total flux in the pho- 
todisk is given by: 

-Fphotodisk = 27rr 2 ,£ . (6) 

The projected surface brightness distribution of the pho- 
todisk inclined at any angle i was calculated by integrating 
Equation 5 over the areas in the face-on photodisk plane 
seen by each projected pixel. The disk component was 
collapsed to a point (zero radius) source anytime its disk 
scale length was less than 0'.'02 (0.2 WFC pixel) when the 
bulge+disk model fits were performed. Equations 3 and 5 
are given in their circularly symmetric form for simplicity, 
but the GIM2D model certainly does not assume circu- 
lar symmetry since it includes photobulge ellipticity and 
photodisk inclination as fitting parameters. 

The optical thickness of disk galaxies remains a hotly de- 
bated issue locally (Valentijn 1994; Giovanelli et al. 1994; 
Disney et al. 1992; Valentijn 1990), and this issue obvi- 
ously has important consequences for the interpretation 
of local and high-redshift photodisk data. The photodisk 
optical thickness is not one of the fitting parameters, but 
photodisks can be fitted with surface brightness profiles of 
different optical thicknesses given by the parameter C a b s 
in the equation: 

Wlphotodisk.obs = ^photodisk,face-on + 2.5 C a bs log(a/6) (7) 

where TO p hotodisk,obs is the observed photodisk total mag- 
nitude and m p hotodisk,facc-on is the face-on photodisk total 
magnitude, a/b is again the axial ratio of the photodisk 
component. Allowable values for C a bs are between (op- 
tically thin photodisks) and 1 (optically thick photodisks). 
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All GSS photodisk models used in this analysis were to- 
tally optically thin (C a bs =0). No internal absorption was 
applied to photobulges. 

A PSF-deconvolved semi-major axis half-light radius 
was also computed for each galaxy by integrating the sum 
of Equations 3 and 5 out to infinity with the best fitting 
structural parameters. This half-light radius may be in 
error for galaxies with large differences between their pho- 
tobulge and photodisk position angles. 

The WFPC2 detector undersampling was taken into ac- 
count by generating the surface brightness model on an 
oversampled grid, convolving it with the appropriate over- 
sampled PSF, shifting its center according to dx and dy 
and rebinning the result to the detector resolution for di- 
rect comparison with the observed galaxy image. Letting 
dx and dy be free parameters in the fits also helps com- 
pensate for possible errors in the initial determinations of 
the galaxy centroid by SExtractor. 

5.3. Types of Fits 

Before describing each type of fit, the question of what 
is being fitted should be discussed. Some studies (e.g., 
Schade ct al. 1996) "symmetrized" galaxy images before 
fitting with a bulge+disk model. Symmetrization removes 
any asymmetric structures around a chosen image pivot 
point, i.e., the galaxy image centroid. Symmetrization 
makes sense since the fitting model is itself symmetric, and 
asymmetries can affect the final values for the best-fitting 
parameters (e.g., Marleau & Simard 1998 and Section 7.3 
of this paper) when these asymmetries are very strong. 
Other studies (e.g., Ratnatunga et al. 1999) do not use im- 
age symmetrization. Although image symmetrization is an 
option in GIM2D, it was not used for the GSS structural 
analysis. To symmetrize or not to symmetrize is really a 
"philosophical" choice. The result of image symmetriza- 
tion depends strongly on the choice of pivot point. An 
error in pinpointing the location of the true pivot point 
may result in an artificial change in the size of the ob- 
ject being fitted. The ambiguity on the true location of 
the galaxy centroid will increase with the strength of the 
asymmetric structure. So, as symmetrization would be- 
come more and more useful, the error introduced by a bad 
pivot point is also prone to become larger. Since sym- 
metrization is not used here, possible effects of asymmet- 
ric structures on fitting parameters should be kept close to 
mind. For example, star-forming regions and spiral arms 
in the disk of galaxy can artificially increase the measured 
disk scale length and decrease the measured bulge fraction 
(Note that the problem can also affect symmetrized images 
if deviant structures were symmetric about the centroid 
of the galaxy). Fortunately, residual image indices (sec- 
tion 5.6) can be used to identify objects for which struc- 
tural parameters may have been seriously compromised 
by asymmetric structures. An improved symmetrization 
procedure in which the total flux removed from the input 
image is minimized might solve the pivot point problem. 
The effects of asymmetries on the measured structural pa- 
rameters from non-symmetrized fits are studied in detail 
in Section 7.3. 

GIM2D offers three types of galaxy surface brightness 
fits: (1) two-bandpass, separate fits, (2) two-bandpass, si- 
multaneous fits, and (3) single-bandpass, multiple image 



fits. The first two types of fits were used for the GSS 
structural analysis, and the third one is particularly useful 
for dithered data such as used for HST/NICMOS. Science 
goals dictate the type of fits that should be used. 

For the separate, two-bandpass fits, the GSS Isu and 
Vqq6 thumbnail images were fitted independently, i. c., no 
fitting parameter was constrained to have the same value 
in both bandpasses. The only connection between the two 
bandpasses was the use of the Ig,u segmentation thumbnail 
image for both Isu and Vqoq fits. By comparing parameter 
values in both bandpasses, this type of fit can provide valu- 
able information about color gradients in photobulges or 
in photodisks. These color gradients would yield very dif- 
ferent photobulge effective radii or photodisk scale lengths 
in the two bandpasses. Also, the difference in the loca- 
tion of the I$i4 and Vena model centroids can be used as 
another image asymmetry estimate. For example, strong 
blue asymmetric structures possibly arising from star for- 
mation should perturb the centroid of the Vqoq model more 
than the I S i4 centroid. The main disadvantage of this type 
of fit is that it docs not make maximum use of all the 
information available at a given signal-to-noise ratio. In 
the absence of color gradients, all the information should 
be used by fitting both bandpasses simultaneously as de- 
scribed below. 

For the simultaneous, two-bandpass fits, the GSS I$i4 
and V(306 thumbnail images were fitted simultaneously with 
all but three fitting parameters forced to take on the same 
values in both bandpasses. The three fitting parameters 
allowed to be different were the total flux in the model 
(now Fj and Fy), the bulge fraction (now (B/T)j and 
(B/T)v) and the background residual level (now dbi and 
dby). This type of fit is, by its nature, blind to color gradi- 
ents in the structural subcomponents, but, in the absence 
of such gradients, it should yield better photobulge and 
photodisk colors. 

The third and last type of fit was not used for the 
DEEP/GSS analysis, but it deserves some discussion for 
its usefulness in fitting stacks of dithered images. Some 
WFPC2 and many NICMOS imaging programs consist of 
dithered exposures of the same region of the sky. The 
offsets between the images can be just a few pixels or 
they may be a significant fraction (~ 1/3) of the detec- 
tor field of view. The offsets are often non-integer pixel 
shifts, and flux interpolation must be used in these cases 
to combine the images. Interpolation will affect the noise 
characteristics of the images and should be avoided. Sig- 
nificant flux errors can also be introduced by interpolating 
undersampled image data (e. g., WF and NICMOS/NIC3 
cameras). The single-bandpass, multiple image fits use the 
same model to simultaneously fit multiple dithered images. 
GIM2D computes an image centroid for each image, and 
the image centroids are passed on to the model generating 
routine so that image models will be shifted accordingly. 
The main advantage here is the ability to avoid interpo- 
lating pixel flux. 

All three types of GIM2D fits are performed on all pixels 
flagged as object or background in the SExtractor segmen- 
tation image. Object areas in the segmentation image are 
sharply delineated by the location of the isophote corre- 
sponding to the detection threshold (1.5-abkg here) since 
SExtractor considers all pixels below this threshold to be 



8 



Simard et al. 



background pixels. However, precious information on the 
outer parts of the galaxy profile may be contained in the 
pixels below that threshold, and fits should therefore not 
be restricted only to object pixels to avoid throwing that 
information away. Pixels belonging to objects in the neigh- 
borhood of the primary object being fitted are masked out 
of the fitting area using the SExtractor segmentation im- 
age. The flux from the primary object that would have 
been in those masked areas in the absence of neighbors 
is nonetheless properly included in the magnitude mea- 
surements given in this paper because magnitudes were 
obtained by integrating the best-fit models over all pixels. 



5.4. Initial image moments and Fitting Parameter 

Values 

Even though the SExtractor local background was sub- 
tracted from each science thumbnail image, GIM2D can 
compute a residual mean background level to correct for 
any error in SExtractor background estimates. GIM2D 
can also be instructed to compute its own estimate of 
o~bkg- To compute background parameters, GIM2D uses 
all the pixels in the science thumbnail image flagged as 
background pixels (flag value of zero) in the SExtractor 
segmentation image. GIM2D further prunes this sample of 
background pixels by excluding any background pixel that 
is closer than five pixels from any (primary or neighboring) 
object pixels. This buffer zone ensures that the flux from 
all SExtracted objects in the image below all the 1.5-abkg 
isophotes does not significantly bias the mean background 
level upwards and artificially inflate Obkg- The GIM2D 
background parameters are tested in Section 7. For the 
GSS fits, background parameters were re-calculated with 
GIM2D before fitting, and the residual background lev- 
els alb were then frozen to their recalculated values in the 
surface brightness fits. 

SExtractor object centroid coordinates can be accepted 
as is, or GIM2D can determine its own intensity-weighted 
object centroid using a multi-threshold, minimum area 
procedure. This procedure first identifies which of three 
thresholds (10-crbkg, 5-o~bk g , and S-abkg) has enough ob- 
ject pixels, as flagged by the SExtractor segmentation im- 
age, to meet a certain minimum area requirement. This 
minimum area was set to 8 for the GSS objects. GIM2D 
then computes intensity-weighted centroid coordinates at 
the highest level with enough pixels. If none of the three 
thresholds provides enough pixels to calculate the centroid, 
then GIM2D simply uses all object pixels to compute the 
centroid. 

It is possible to place generous initial limits on certain 
structural parameters (total flux, initial photobulge and 
photodisk scale radii and position angles) based on simple 
image moments widely in use throughout the literature. 
For a given object, these moments are computed about 
the image centroid over the object pixels. These moments 
include: 
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where the sums are over the total number of pixels belong- 
ing to the object, Xij and yij are pixel coordinates with 
respect to the object centroid, and is the background- 
subtracted pixel flux in the (i, j)th pixel. The initial total 
flux estimate for the image model is given by Equation 8a, 
and the maximum limit on this total flux is set to twice 
M to t- The initial values for the photobulge effective radius 
and photodisk scale length are both set to the intensity- 
weighted average radius of Equation 8e, and the maximum 
limits placed on these two scale radii are set to twice the 
intensity- weighted average radius M rr . The minimum lim- 
its on the two scale radii are set to zero to allow the model 
to fit unresolved sources if needed. 

The object position angle is another parameter that can 
be set to an initial value given by the image moments above 
as: 
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(9) 



So, both 4>b and (f>d are initially set to <f>. No mini- 
mum/maximum limits arc placed on the photobulge and 
photodisk position angles since Equation 9 is a circu- 
lar function. In the case of a perfectly circular object 
(M xx — M yy ), both <f>b and (f>d are set to zero. Both param- 
eters are given generous initial Metropolis "temperatures" 
(see next section) of 60°. 

5.5. Metropolis Fitting Algorithm 

The 12-dimensional parameter space can have a very 
complicated topology with local minima at low S/N ra- 
tios. It was therefore important to choose an algorithm 
which did not easily get fooled by these local minima. 
The Metropolis algorithm (Metropolis et al. 1953; Saha & 
Williams 1994) was designed to search for optimal param- 
eter values in a complicated topology, and it is widely used 
in many fields of physics and computer science. Compared 
to gradient search methods, the Metropolis algorithm is 
not efficient, i. e., it is CPU intensive. On the other hand, 
gradient searches are "greedy." They will start from ini- 
tial parameter values, dive in the first minimum they en- 
counter and claim it is the global one. 

The first step in the fitting process is the "initial con- 
dition filter" (ICF). The ICF coarsely samples the very 
large volume of structural parameter space Vo defined by 
broad limits from the initial image moments (Section 5.4) 
and user-specified parameter constraints to determine a 
promising sub-volume to be explored by the Metropolis 
algorithm. The uniform ICF sampling for the i th struc- 
tural parameter follows the equation: 



Xi,o + (u - \)Ti 



(10) 

where Axi is a trial step in parameter space, x^o is the 
i th coordinate of the ICF sampling origin given by the 
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image moments and user constraints, u is a uniform de- 
viate between and 1, and Tj is the "Metropolis temper- 
ature" of the parameter. The Metropolis temperature of 
a parameter is given in the same units as that parame- 
ter. The ICF samples parameter space with Njcf (set to 
400 here) different galaxy image models and keeps track 
of which model yields the highest likelihood. The appro- 
priate value of Njcf depends on how well constrained the 
search volume is initially. After completion, the ICF sets 
the sampling origin to the parameters of that best ICF 
model. The volume Vb is reduced by a factor of Njcf by 
cooling each Metropolis parameter temperature according 
to T! = T l /{N IC F) 1/n with n = 12 being the number of fit- 
ting parameters in the galaxy image model used here. The 
new sampling origin, the reduced search volume, and the 
cooler Metropolis parameter temperatures are then sent to 
the Metropolis algorithm to continue the search further. 

The Metropolis algorithm in GIM2D starts from the 
sampling origin given by the ICF and computes the like- 
lihood P(w\D, M) that the parameter set w is the true 
one given the data D and the model M. It then generates 
random perturbations Ax about that initial location with 
a given "temperature." When the search is "hot," large 
perturbations are tried. After each trial perturbation, the 
Metropolis algorithm computes the likelihood value Pi at 
the new location, and immediately accepts the trial per- 
turbation if Pi is greater than the old value Pq. However, 
if Pi < P , then the Metropolis algorithm will accept the 
trial perturbation only P\/Pq of the time. Therefore, the 
Metropolis algorithm will sometimes accept trial perturba- 
tions which take it to regions of lower likelihood. This ap- 
parently strange behavior is very valuable: if the Metropo- 
lis algorithm finds a minimum, it will try to get out of it, 
but it will only have a finite probability (related to the 
depth of the minimum) of succeeding. The "temperature" 
is regulated according to the number of accepted itera- 
tions. If the Metropolis accepts too many trial perturba- 
tions, then the search is too "cold," and the temperature 
must be increased. Conversely, if the Metropolis rejects 
too many trial perturbations, then the search is too "hot," 
and the temperature must be decreased. The Metropolis 
temperature is regulated such that half of the trial pertur- 
bations are accepted. The temperature is checked every 
fiftieth iteration, and the terms of the covariance matrix 
s are adjusted according to whether the search is too hot 
or too cold. The more commonly known "simulated an- 
nealing" technique is a variant and a special case of the 
Metropolis algorithm in which the temperature is only al- 
lowed to decrease until a "ground-state" is reached. 

The step matrix for the trial perturbations Ax is given 
by the simple equation Ax = Q • u, where the vector u 
consists of randomly generated numbers between —1 and 
1. The matrix Q thus controls the step distribution, and 
random steps with any desired covariance can be gener- 
ated by solving the equation s = Q • Q through Cholcski 
inversion. The matrix s is the local covariance matrix of 
accepted iterations (Vanderbilt & Louie (1984)). In short, 
the Metropolis sampling of parameter space shapes itself 
to the local topology. 

Convergence is achieved when the difference between 
two likelihood values separated by 100 iterations is less 
than 3er of the likelihood fluctuations. After convergence, 



the Metropolis algorithm Monte-Carlo samples the region 
where the likelihood is thus maximized and stores the ac- 
cepted parameter sets as it goes along to build the distri- 
bution P(w\D, M). Once the region has been sufficiently 
sampled (N samp i e = 300 here), the Metropolis algorithm 
computes the median of P(w\D,M) for each model pa- 
rameter as well as the 68% confidence limits. The output 
of the fitting process consists of a PSF-convolved model 
image O, a residual image R and a log file containing all 
Metropolis algorithm iterations, the final parameter values 
and their confidence intervals. 

GIM2D creates two output images for each fitted object: 
an image of the PSF-convolved model and an image show- 
ing the residuals from the bulge+disk model subtraction. 
Figure 3 shows mosaics (science, GIM2D output models, 
and residual) of 35 GSS galaxies with igi4 < 24 chosen 
at random from the whole GSS sample. All three mosaics 
use the same greyscale levels. The Igu magnitude, bulge 
fraction and semi-major axis half-light radius are given on 
the original images. 

5.6. Residual Image and Asymmetry Indices 

GIM2D computes six image indices from the thumbnail 
residual image that can be used to globally characterize 
the structures left after the best galaxy image model has 
been subtracted. 

The first indices, Rt and Ra, were first applied to dis- 
tant galaxies by Schade et al. (1995). These indices are 
based on local studies of spiral arm patterns (Elmegreen, 
Elmcgrccn, & Montenegro 1992), and they provide an esti- 
mate of the overall smoothness of the galaxy image with re- 
spect to the fitting model. Following Schade et al. (1995), 
Rt and Ra are defined as: 
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where the Rij's are pixel values in the residual image, and 
the i?^ 80 's are pixel values in the residual image rotated by 
180°. Since these raw values (RT)raw and (Ra)™™ involve 
taking absolute values of pixel fluxes, they will yield a pos- 
itive signal even in the sole presence of white noise. These 
raw values must therefore be background-corrected. The 
background corrections, (Rr)bkg and (RA)bk g , are com- 
puted over pixels flagged as background pixels in the SEx- 
tractor segmentation image. The Bij's are background 
pixel values in the residual image, and the Bjj 0, s are back- 
ground pixel values in the residual image rotated by 180°. 
The background corrections are computed over a back- 
ground pixel area equal to the pixel areas over which the 
raw indices were computed. Bij's were randomly drawn 
from the full pool of all background pixels in the thumbnail 
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science image to decrease the vulnerability of background 
parameter determination to possible localized background 
image artifacts. 

There is a fundamental difference between the Rt and 
Ra indices of Schade et al. (1995) and the ones imple- 
mented in GIM2D. Schade et al. calculated residual in- 
dices within a physical radius of 5 kpc irrespective of the 
physical size of the whole galaxy. On the other hand, 
GIM2D computes its Rt and Ra indices within ten circu- 
lar apertures whose radii are multiples of the galaxy half- 
light radius ru ranging from 1 ru to 10 rta- The GIM2D 
Rt and Ra indices therefore sample the same fractional 
area for all galaxies. Schade et al. define normal galaxies 
as galaxies with Rt + Ra < 0.14. Im et al. (2001) were 
able to define a sample of high-redshift E/SO purely based 
on quantitative morphology with only two simple selection 
criteria ((-B/T)/ > 0.4 and R T + R A < 0.08), and Mcin- 
tosh (2001) used Rt + Ra to study the SO populations of 
a sample of local Abell clusters. Measuring Rt and Ra 
in different bandpasses could also show the wavelength- 
dependence of the strength of residual structures. If Rt 
and Ra have larger values in bluer bandpasses (e. g., V<306 
versus Isu), then this would suggest that asymmetries 
arise from star-forming regions. For example, there ap- 
pears to be a correlation between the strength of residual 
structures and star formation rates measured from [Oil] 
emission lines in galaxies in poor groups (Tran et al. 2001). 

The next two image indices are from the automated clas- 
sification proposed by Abraham et al. (1994, 1996). This 
classification relics on two parameters: one measuring the 
concentration of galaxy light (C) and the other one mea- 
suring image asymmetry (A). The so-called C — A sys- 
tem explicitly takes the ellipticity of galaxy images into 
account instead of simply using circular apertures. The 
image moments are used to define an equivalent elliptical 
distribution. For each object, the area Ai a of the 2-abkg 
isophote is first computed. Then, following Abraham et 
al. (1994), the definitions of the concentration index C 
and normalized radius rij are: 
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where E(rtj) is an elliptical aperture bounded by ry, i. e., 
is constant on elliptical isophotes, and the image mo- 
ments are normalized such E(l) = A^ a - Abraham et al. 
adopted a value of a = 0.3. C(a) was computed for four 
values of a (0.1, 0.2, 0.3, and 0.4) for the GSS galaxies. 
The A parameter measures image asymmetry in a similar 
way to Ra- However, A is computed directly on the science 
image whereas Ra is computed on the residual image from 
the bulge+disk fit. A is also computed over the area A 2cr , 
and again similarly to Ra, a background correction must 
be applied to raw A values to remove the extra positive 
signal from the background image noise. The C — A indices 
have been implemented in GIM2D exactly as prescribed by 
Abraham et al. (1996). Subsequent works (e.g., Conselice, 
Bershady, & Jangren 2000; Wu 1999) have since shown 



that choosing an image pivot point which minimizes the 
measured asymmetry greatly improves the classification 
results, but these improvements have not been included in 
GIM2D yet. 

In addition to the previous indices, two new indices, 
A z and D Zl were defined to quantify residual asymmet- 
ric structures. The A z index is calculated over the pixels 
belonging to the object or the background. Each pixel is 
compared to its symmetrical counterpart about the center 
of the object. If the flux in that pixel is no~bkg higher than 
the flux in its symmetrical counterpart, then the pixel is 
taken to be part of an asymmetrical component. A z is the 
sum of the fluxes of all such pixels normalized by the total 
object flux. The A z index is computed for n = 2, 3 and 5 
within circular annuli with radii between lr^i and lOr^/. A 
statistical background correction similar to the ones used 
for (Rt) and (Ra) was applied to the raw A z values. The 
D z index takes advantage of the isophotal shape of the 
object as measured by SExtractor. D z is the sum of the 
fluxes of the object pixels (as given by the segmentation 
image) with symmetrical counterparts about the object 
center which do not belong to the object. This sum is nor- 
malized by the total object flux. D z is computed outside 
of one half-light radius, and it is sensitive to asymmetries 
such as tidal arms. 

6. STRUCTURAL CATALOGS 

The images of GSS galaxies were fitted with both the 
separate and simultaneous fitting procedures, and exten- 
sive sets of GSS simulations were performed for both (Sec- 
tion 7). As a result, four structural parameter catalogs are 
presented in this paper: two science catalogs and two sim- 
ulation catalogs. These catalogs are automatically gener- 
ated from the GIM2D output log files by a set of scripts. 
These scripts calculate all the final parameter values in- 
cluding physical lengths and rest-frame quantities (Sec- 
tion 6.2 below) with full Monte-Carlo propagation of the 
parameter errors (Section 6.3 below). The Keck/LRIS 
spectroscopic redshifts used to compute physical radii in 
kiloparsecs and rest-frame magnitudes and colors for the 
whole galaxy, the photobulge and the photodisk are taken 
from Phillips et al. (2002). The observational, reduction 
and analysis procedures used to measure these redshifts 
are fully described in their paper. 

6.1. Description 

The contents of the science catalogs are listed in 
machine-readable Tables ?? (separate fits) and ?? (simul- 
taneous fits). The separate fit science catalog contains 
7450 objects, and each object has 330 columns of parame- 
ter information (including error bar columns). Simultane- 
ous fits were performed only to objects with Keck/LRIS 
redshifts from the DEEP survey of the Groth Strip. The 
simultaneous catalog thus contains 648 objects with 279 
columns of information for each one, including error bar 
columns. There are two error columns associated with 
each entry described in Tables ?? and ??: one for the 
lower 68% confidence bound and the other for the upper 
68% bound. Most of the column descriptions in Tables ?? 
and ?? are self-explanatory, but some of them require 
further details: 
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DEEP '/GSS IDs "gssid": The internal DEEP/GSS ob- 
ject IDs (" gssid') are given by the format FFC-XXYY 
where FF is the GSS field (Table 1, Column 1), C is the 
WFPC2 chip number , and XX and YY are the object co- 
ordinates on the chips in units of 10 pixels. These internal 
IDs are extended with letters ("a", "b", "c", etc.) when a 
group of objects are close enough together that they would 
all have the same primary ID. In addition to these inter- 
nal DEEP/GSS objects ID's, the catalogs list the J2000.0 
coordinates of each object. 

Physical scale lengths: All angular scale lengths (r^i, r e , 
and Td) were converted to physical lengths according to 
the equation: 



lOOOffol + ^o v/f! m (l + z') 3 + 
where z is the redshift, c is the speed of light, r is the 
measured angular scale length in radians, and R is the 
corresponding physical length in kiloparsecs (Hogg 1999). 
Equation 13 is only valid for flat (Clk = 1 — Q m — = 0) 
cosmologies. 

Position angles on sky: The GIM2D photobulge and 
photodisk position angles 4>b,image and 4>d,image measured 
on the images clockwise with respect to the positive y-axis 
were converted to real position angles <fib,sky and <j>d,sk y on 
the sky using the telescope position angles stored in the 
HST image headers. 

AB Magnitudes: Even though Vega-based magnitudes 
are primarily used in this paper, Tables ?? and ?? also pro- 
vide AB magnitudes. To go from Vega to AB magnitudes, 
one uses 7 8 i 4 (AB) = 7 8 i 4 (Vega) + 0.44 and V 6 06(AB) = 
V 606 (Vega)+0.11. 

Galaxy rest-frame B -band photobulge fraction: The rest- 
frame photobulge fraction is given by the simple equation: 

(B/T) Barest = 10.0 (MB - galaxy - MB - bulBo)/2 ' 5 (14) 

where Ms jga i axy and M^buige are the rest-frame S-band 
magnitudes of the galaxy and photobulge respectively. 
Since different fc-corrections apply to the photobulge and 
photodisk stellar populations, the observed bulge fraction 
of a galaxy will change with redshift. So, rest-frame pho- 
tobulge fractions make more uniform photobulge fraction 
selections possible. 

6.2. Rest-Frame Quantities 

Rest-frame quantities in the GIM2D/GSS structural 
catalogs were calculated using two very different sets of fc- 
corrections. These two sets provide independent checks of 
the reliability of the resulting rest-frame quantities. Rest- 
frame quantities were calculated independently for the to- 
tal galaxy, the photobulge and the photodisk. 

The first set of fc-corrections (referred to as Gronwall fc- 
corrections throughout this paper) is based on the work 
of Gronwall & Koo (1995). The corrections are based 
on eleven theoretical galaxy spectral energy distributions 
from the 1995 Bruzual and Chariot models. The param- 
eters of these theoretical SEDs are given in Gronwall & 
Koo (1995). Some of these SEDs include dust and star- 
bursting populations. The rest-frame magnitudes and col- 
ors of any galaxy are obtained by interpolating the SEDs. 
The input quantities are galaxy redshift, Vega-based I$i4 



magnitude and V<306 — -^814 color. The output consists of 
rest- frame BVRI absolute magnitudes, rest-frame (U — B) 
and (B — V) colors, and rest-frame B and K magnitudes. 

The second set of fc-corrections (referred to as Willmer- 
Gcbhardt or WG fc-corrections throughout this paper, 
Gcbhardt et al. 2002) is based on actual galaxy spectra. 
These spectra are taken from the Database of UV-Optical 
spectra of Nearby Quiescent and Active Galaxies (Kinney 
et al. 1996; Schmitt et al. 1997). This database has re- 
cently been expanded to include 99 galaxies, 48 of which 
have full wavelength coverage from 1200 to 10000A with a 
combination of International Ultraviolet Explorer (IUE) 
and ground-based spectra. Filter bandpasses are con- 
volved with the galaxy spectra to produce rest-frame 
(U — B) colors and -B-band fc-corrections k B - A polyno- 
mial is fitted to the observed (V<306 — hu) colors and the 
rest-frame (U — B) SED colors in each redshift range. The 
best-fit polynomial reproducing the data over the redshift 
range 0.1 — 1.1 is given by: 

(U - B) WG = -0.8079 - 0.049752z - 1.6232z 2 + 1.04067z 3 
+ 1.5294z 4 - 0.41190x 5 - 0.56986z 6 + (0.61591 
+1.072492! - 2.2925z 2 + 1.3370z 3 )(V 606 - J 8 i 4 ) 
+ (0.280481 - 0.387205z + 0.043121z 2 ) 
(V 606 - Isu) 2 (15) 
Similarly, the B-band fc-corrections are given by: 

fc B j = 0.0496 + 0.46057z + 1.40430z 2 - 0.19436z 3 

-0.2232z 4 - 0.36506z 5 + 0.17594z 6 + (2.0532 
-2.8326z + 1.05580z 2 - 0.67625z 3 )(V 606 - hu) 
+ (0.10826 - 0.68097z + 0.61781z 2 )(V 606 - huf, 

(16) 

and the rest-frame B-band magnitude in the Willmer- 
Gebhardt system is given by: 

m b ,wg = hn - DM(n m ,n A ,n k ) + k BI (17) 

where DM(Q m , Oa, ^k) is the distance modulus for the 
adopted cosmology. See Gebhardt et al. (2002) for more 
details. 

6.3. Error Estimates 

All parameter error estimates in the GIM2D/GSS struc- 
tural parameter catalogs are 68% confidence limits. Many 
of these error estimates are asymmetric since they were 
derived through full Monte-Carlo propagations of the pa- 
rameter probability distributions P(w\D, M) computed 
by GIM2D through all the transformations required to 
calculate a given final parameter. This process takes into 
account all the Gaussian and non-Gaussian covariances 
among the parameters. To illustrate the process, consider 
the observed photobulge I$u magnitude of a galaxy. This 
quantity depends on both the total galaxy model flux and 
the observed F81AW bulge fraction. So, the parameter 
probability distributions P(F 814 \D, M) and P(B/T\D, M) 
were first Monte-Carlo sampled 500 times, and a photo- 
bulge magnitude was calculated each time using Equa- 
tion 1. The resulting 500 photobulge magnitudes were 
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then sorted, and the lower and upper 68% confidence er- 
ror estimates were derived from that sorted distribution. 
Going one step further, the Vqo6 photobulge magnitude 
can be computed in exactly the same way as Isi4, and the 
resulting 500 Veoe ~ -^814 colors can then be transformed 
to photobulge rest-frame B-band magnitudes using Equa- 
tions 16 and 17 so that one can in turn compute the 68% 
lower and upper error estimates on the photobulge abso- 
lute magnitude. 

7. SIMULATIONS 

Three sets of GIM2D simulations were run to char- 
acterize the systematic biases and random errors in the 
GIM2D/GSS structural measurements. These simulations 
are a key element in the interpretation of the observations, 
and the simulation catalogs are presented in the same way 
as real catalogs to emphasize this point. Given any real sci- 
ence plot, it is straightforward to select sets of simulated 
galaxies with the same selection criteria as the observa- 
tions to immediately evaluate the biases and errors present 
in the data shown on that plot. The first set of simulations 
applies to separate fits of 7si4 and Vqo6 GSS galaxy images 
(Section 7.1). It contains 5995 simulations, comparable in 
size to the real science catalog. The second set includes 
5195 GSS simultaneous Vqoq/Isu fit simulations, and these 
simulations cover the full range of observed GSS photob- 
ulge and photodisk V<306 — ^814 colors (Section 7.2). The 
two sets of simulations above only include smooth galaxy 
image simulations. The effects of asymmetric structures on 
the measured structural parameters were explored with a 
third set of simulations described in Section 7.3. 

7.1. Separate Fits 

For the GSS separate fit simulations, 5995 smooth 
galaxy image models were created with structural pa- 
rameters uniformly generated at random in the follow- 
ing ranges: 20.0 < I 814 < 25.0, 0.0 < B/T < 1.0, 
< r e < 0'.'7, 0.0 < e < 0.7, < r d < 0'.'7, and 
< i < 85°. The Sersic photobulge index was held fixed at 
n = 4 for all models. Both photobulge and photodisk po- 
sition angles were fixed to 90° for all simulations, and the 
bulge and disk sizes were uniformly generated in the log 
of the size ranges above. The goal of these simulations is 
to characterize biases and errors and not to simulate what 
the real Universe would look like through the GIM2D ob- 
servational "lens." The uniformity of the parameter dis- 
tributions adopted here is therefore perfectly suitable to 
the task even though real galaxy parameters (e. g., bulge 
fraction) may not be so distributed. In the same spirit, no 
correlations were imposed between the input parameters 
despite the fact that some parameters (e. g., r e and rd) 
may be correlated in some types of galaxies (e.g., Courteau 
et al. 1996). 

Each simulation was convolved with a F814W TinyTim 
PSF. This PSF had the same parameters as the TinyTim 
F8UW PSFs used in the GSS analysis (Section 5.1)! The 
same PSF was used in both creating and analyzing the 
simulations, so the results will not include any error in the 
structural parameters due to PSF mismatch. Poisson devi- 
ates were used to add photon noise due to galaxy flux into 
the simulations. The noisy images were then embedded in 
a 20" x 20" section of one of the real F81AW GSS images 



to provide a real background for the simulations. In addi- 
tion to sky photon noise and detector read-out noise, the 
real background noise includes brightness fluctuations of 
very faint galaxies below the detection threshold. The sim- 
ulations were SExtracted with exactly the same SExtrac- 
tor parameter files (Sections 4.1 and 4.2) as used for the 
GSS analysis, and GIM2D extracted science and segmen- 
tation thumbnails from the simulations following exactly 
the same steps as for the real galaxies (Section 4.3). Fi- 
nally, the GIM2D output log files were processed through 
the same scripts to produce a catalog of final recovered 
structural parameters. The content of this catalog is listed 
in Table ??. 

7.1.1. Systematic and Random Errors 

For the sake of simplicity, the main tool adopted here 
to visualize errors is a set of two-dimensional maps giving 
systematic and random errors at each position. It should 
therefore be kept in mind that these maps can only offer 
a limited representation of the complex multidimensional 
error functions. As the large number of parameters in Ta- 
bles ?? and ?? indicates, a full description of all system- 
atic and random errors over all of bulge+disk multivariate 
structural space would considerably add to the length of 
this paper. Therefore, the error analysis presented in this 
section will focus on only three main galaxy structural pa- 
rameters: total apparent magnitude, bulge fraction and 
half-light radius. Errors on any other set of parameters 
can be described in the same way, and the simulation cat- 
alog can be used to tailor error analyses to the needs of 
the specific science goals being pursued. 

The error maps can be cast in terms of input or re- 
covered coordinates, and the choice of coordinate system 
depends on how the error maps will be used. Input coor- 
dinates (i. e., the "true" coordinates) can be used to com- 
pute errors that are to be applied to theoretical galaxy 
structural catalogs in order to convert them to observed 
quantities (e.g., Simard et al. 2002. To illustrate this pro- 
cess, let wt be the position of a mock galaxy in theoretical 
structural space, and let rhi,T be its theoretical half-light 
radius. If the simulation catalog shows that the recovered 
half-light radii of galaxies at wf are systematically in error 
by an amount Ari, then let r' hl T = r^^r + Ar x . This new 
radius r' hl T is not yet the same as an observed radius as 
it does not include a random error. The random error on 
the half-light radius o(rhi,T) at wf can also be calculated 
from the simulation catalog, and another radius correc- 
tion Ar2 drawn at random from a Gaussian distribution 
of width cr(rhLT) can be applied to r' hl T to produce the 
final "observed" theoretical half-light radius. The second 
set of coordinates, the recovered quantities, is the simula- 
tion equivalent of observed quantities, and error maps cast 
in those coordinates can be directly compared to the real 
data to see how important errors are in different regions 
of the observational space. 

Figures 4, 5, and 6 show maps of errors on the galaxy to- 
tal magnitude I^u, galaxy half-light radius rhi and galaxy 
bulge fraction (B/T) respectively as a function of galaxy 
magnitude and galaxy half-light radius for the DEEP/GSS 
separate structural fits. The two top panels in each fig- 
ure show the mean parameter error (left-hand panels, 
top number in cells) and the \-a parameter random er- 
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ror (right-hand panels, top number in cells) as a function 
of input galaxy magnitude and size. Each cell also gives 
the number of simulations created for that cell (bottom 
number). The simulations are not evenly distributed over 
the galaxy magnitude-log size plane since the simulations 
were uniformly generated in log r e and in log r^. The 
lower left-hand panels show the mean parameter error as 
a function of recovered galaxy magnitude and size, and 
the lower right-hand panels of the three figures show the 
actual DEEP/GSS magnitude-size data. The unresolved 
objects (log Thi,obs < —1-5) are nicely separated from the 
galaxies in these panels. 

Figure 4 shows that the systematic errors on Jgi4 galaxy 
magnitudes start to become significant (A/§i4 ~ —0.2) 
fainter than 7gi4 — 23.5, and that, at a given magnitude, 
errors are larger for the largest galaxies in the simulations 
(log rhi^input > —0.25). The random magnitude errors are 
about 0.05 for I sl4 < 23.0 and 0.13 for I S u > 23.0. The 
systematic errors cast in terms of input or recovered co- 
ordinates are essentially the same since the errors are not 
large enough compared to the cell sizes (0.5 mag and 0.5 
log rhi) to shift simulations from cell to cell. If galaxies 
were pure bulges or pure disks, then according to equa- 
tions 4 or 6, systematic magnitude and size errors should 
be anti-correlated, i. e., given an observed surface bright- 
ness profile and an underestimate, say, of the total flux 
(positive magnitude error), the profile modelling should 
try to compensate for that magnitude underestimate by in- 
troducing a negative error in the size. This anti-correlation 
should still hold for composite systems. A comparison of 
Figures 4 and 5 does show that magnitude and size errors 
are typically anti-correlated. Magnitude and size errors 
are important for fainter and larger galaxies since their 
relatively low surface brightness makes them more vulner- 
able to sky estimate errors. 

The error maps for the galaxy bulge fraction (Figure 6) 
show that bulge fractions are underestimated by about 
0.15 at magnitudes fainter 7si4 = 23.5 with random er- 
rors around 0.25. However, Figure 6 is not really the 
best way to truly understand the behavior of the recovered 
bulge fractions. There are in fact two expected biases in 
the bulge fractions, and these biases arise from two ingre- 
dients of the DEEP/GSS bulge+disk analysis: (1) bulge 
fractions are constrained to stay between and 1, and (2) 
bulge+disk models were fitted to all detected objects irre- 
spective of the signal-to-noise ratio (S/N) of their images. 
The constraint on the bulge fraction forces the recovered 
bulge fractions of both very low (B/T ~ 0) and very high 
(B/T ~ 1) systems to scatter above zero and below one, 
and this bias will affect all galaxies irrespective of their 
S/N ratios. The second bias is inherent to bulge+disk 
model fits to objects with different S/N ratio. Previous 
studies have adopted a two-tier approach to this problem. 
Schade et al. (1995, 1996) first fit pure bulge or pure disks 
to their objects and then decide upon visual inspection of 
the residuals whether a bulge+disk model would be more 
appropriate. Ratnatunga et al. 1999 fit bulge+disk mod- 
els to objects above a certain signal-to-noise, and objects 
below that threshold are fitted only with either a pure 
bulge or a pure disk model. 

A different approach was taken here to deal with the 
bulge+disk S/N bias. Bulge+disk models were fitted to 
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all detected objects here for simplicity and for the sake of 
producing homogeneous structural catalogs. Bulge+disk 
models will converge to a pure bulge or a pure disk model 
only when the signal-to-noise ratio is high enough to defi- 
nitely establish the presence of one and only one structural 
component. At low S/N ratios, the model will always be 
able to "slip" in both structural components. For example, 
the model could make use of a very large disk component 
to compensate for an underestimate of the sky level that 
may have been computed during the fit to a pure bulge 
system, and this would artificially decrease the recovered 
bulge fraction. Low signal-to-noise can also be responsi- 
ble for hiding the outer wing of steep surface brightness 
profiles such as the r 1 / 4 profile into the background noise 
and thus making them harder to identify. Figure 7 best 
shows the bulge fraction biases for the DEEP/GSS sepa- 
rate structural fits. Figure 7 is similar to previous error 
maps except that the errors are now computed over input 
bulge fraction and input galaxy magnitude instead of mag- 
nitude and size. As expected, the B/T systematic errors 
in Figure 7 show that (1) B/T is indeed overestimated 
in the first B/T bin, (2) B/T is underestimated in the 
last B/T bin, and (3) the magnitudes of the discrepan- 
cies increase with magnitude. The homogeneous approach 
to bulge+disk model fitting adopted here is valid as long 
as the results are used in conjunction with careful error 
characterization from the simulation catalogs. 

7.2. Simultaneous Fits 

The GSS simultaneous fit simulations use the isi4 sim- 
ulations of Section 7.1 as a starting point. A compan- 
ion V(506 simulation was created for each I 8 i4 image with 
the same structural parameters except for total flux and 
bulge fraction. The Vqoq total flux and bulge fraction 
were calculated from the F814W total flux and bulge frac- 
tion and from randomly generated photobulge and pho- 
todisk (Ve;o6 — -^8i4) colors. The photobulge and photodisk 
(V(so6 — -^814) colors were uniformly and independently gen- 
erated in the range 0.5-2.2. This range of colors spans the 
full range of observed colors out to a redshift of z = 1.1 in 
the DEEP/GSS survey (Phillips et al. 2002), and it allows 
one to study the effects on fitting results of differences in 
photobulge and photodisk colors. For Jgi4 pure photob- 
ulge systems ((B/T)/ = 1.0), the F606W bulge fraction 
and total flux are given by: 

(B/T) v = 1.0 (18a) 

F tot v = F tot / ^l (l-26-(V-/) photobulgo )/2.5 (18b) 
tl 

where F totJ and F toty are total F814VF and F606W 
galaxy model fluxes in DU respectively, t\ and ty are 
the F 8UW and F606W total exposure times (4400 sec- 
onds and 2800 seconds), and (B/T)/ and (B/T) v are the 
F&1AW and FQ06W photobulge fractions. The zeropoint 
difference between Veoe and 7si4 is 1.26. 

For 7 8 i4 pure photodisk systems ((B/T)j = 0.0), the 
F606W bulge fraction and total flux are given by: 

(B/T) v = 0.0 (19a) 
F toty = FtotJ hL 1Q (^-(v-i) phatad ^)/2. 5 (19b) 
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For Jgi4 composite galaxy systems (0 < (B/T)j < 1), 
the FQ06W bulge fraction and total flux are given by the 
equations: 

{B/T)v = ( 1 7^r loA(y ~ /)/2 ' 5 + L0 ) 1 (20a) 

F tot , v = F toM g^^lO( 1 - 26 -^^)— )/2.5 

(20b) 

where A(V - I) = (V - I) p hotobui g o - (V - I) p hotodisk- 

After adding in Poisson noise, the V<306 simulations were 
also embedded in the corresponding 20" x 20" section of 
one of the real F606T4 7 GSS images. This section of the sky 
was identical to the one used for the /§i4 simulations. As 
was done for the separate fit simulations, the simultane- 
ous fit simulations were processed in exactly the same way 
as the real galaxies to produce a catalog whose content is 
listed in Table ??. Note that, as for the observations, the 
I$i4 segmentation thumbnail images were used for both 
bandpasses in the simultaneous fits. 

7.2.1. Systematic and Random Errors 

Figures 8, 9, and 10 show maps of errors on the galaxy 
total magnitude Isu, galaxy half-light radius r^i and 
galaxy bulge fraction (B/T) respectively as a function 
of galaxy magnitude and galaxy half-light radius for the 
DEEP/GSS simultaneous structural fits (see Figures 4, 5, 
and 6 for the corresponding separate fit results). The two 
top panels show the mean parameter error (left-hand pan- 
els, top number in cells) and the l-cr parameter random er- 
ror (right-hand panels, top number in cells) as a function of 
input galaxy magnitude and size. Each cell also gives the 
number of simulations created for that cell (bottom num- 
ber). Again, the simulations are not evenly distributed 
over the galaxy magnitude-log size plane since the simu- 
lations were uniformly generated in log r e and in log r^. 
The lower left-hand panels show the mean parameter error 
as a function of recovered galaxy magnitude and size. The 
lower right-hand panels of the three figures show the ac- 
tual DEEP/GSS magnitude-size data. The number of data 
points is not nearly as large as in Figures 4, 5, and 6 since 
simultaneous bandpass structural fits were performed only 
on DEEP/GSS galaxies with secure Keck/LRIS rcdshifts. 
The errors from the simultaneous structural fits behave 
the same way as the errors from the separate fits, and si- 
multaneous fit errors seem to be slightly smaller than those 
from the separate fits as one would expect from simultane- 
ously using all the information content of both bandpasses. 
However, the improvement in the errors may not be as 
marked as expected since the simultaneous fit simulations 
included varying bulge fraction and colors as additional 
input parameters. 

Bulge and disk colors were included in the simultane- 
ous structural fit simulations with the goal of testing how 
well they are recovered in the fits. Figure 11 shows the 
structural component colors recovered by GIM2D from 
simultaneous structural fit simulations for galaxies with 
rw,8i4 > 0"15, and r hl ,606 > 0'T5. The V 606 limits for the 
galaxy, bulge, and disk magnitudes were set to 26.0. The 



nine panels show recovered Vqoq — Isu colors versus input 
^606 — -^814 colors for the galaxy as a whole (top panels), 
the bulge (middle panels) and the disk (bottom panels) 
in three different magnitude ranges. The colors are all 
well recovered by the fits. The mean and rms color dif- 
ference in the three magnitude ranges are (0.002, 0.020), 
(0.010, 0.034), and (0.063, 0.087) for the galaxy as a whole, 
(-0.001, 0.255), (0.022, 0.120), and (0.034, 0.186) for the 
bulge and (0.008, 0.047), (0.006, 0.211), and (0.030, 0.304) 
for the disk. There are no significant systematic color er- 
rors, and the rms scatter increases with magnitude. 

Although the recovered colors show no systematic off- 
sets from the input colors, there are interesting outliers 
in some panels of Figure 11. In the leftmost middle 
panel, some recovered bulge colors are much too blue 
compared to their input colors. The three most dis- 
crepant bulges ((V 606 - hiijinput > 1-9 and (V e0 e - 
1 814 ) recovered < 1-3) are all very red bulges with very blue 
disks {{V 606 - hu) input, disk < 0.90), and their effective 
radius differs from the scale length of their disk by a fac- 
tor of five or more. In contrast, the central middle panel 
(21 < Isu{bulge) < 22.0) shows recovered bulge colors 
which are too red for their input colors. The two outlier 
bulges ((V606 — hu)in P ut < 1-0) with red recovered col- 
ors ((V606 — Isu) recovered > 1-5) are quite blue compared 
to their disks ((V 606 - 4u) input, disk > 2.1). One of the 
bulges has an effective radius that differs by a factor of 
10 from the disk scale length, but the other has an ef- 
fective radius comparable to the disk scale length. The 
central bottom panel shows recovered disk colors that are 
too blue compared to their input values. The two right- 
most bottom outliers are very red disks with bluer bulges 
((Veoe - hu) input, bulge < 1-0), and their effective radius is 
different from the disk scale length by a factor of 20-23! 

The outliers in Figure 1 1 lead to an important question: 
Is there a combination of bulge fraction, bulge/disk size 
ratio r e /rd, and bulge/disk colors for which bulges can 
be mistaken for disks and vice versa? Figure 12 shows 
the systematic error (mean error) on bulge fraction as a 
function of input bulge fraction and input log bulge/disk 
size ratio (r e /r<j) for both the separate structural fit (SPF) 
simulations and the simultaneous structural fit (SMF) sim- 
ulations. There are no regions of that plane in which bulge 
fractions are systematically in error. This confirms the ab- 
sence of systematic deviations in bulge and disk colors in 
Figure 11. However, there are places (especially for very 
small or very large bulge/disk size ratios) where the mini- 
mum or the maximum bulge difference is quite large, and 
these extrema can account for the kind of color outliers 
seen in Figure 11. 

7.3. Effects of Asymmetric Structures On Fitting 
Parameters 

Non-smooth local features in a galaxy 2D light profile 
can alter the best parameters derived with GIM2D de- 
pending on their brightnesses and positions in the galaxy. 
For example, a very bright feature at the center of the 
galaxy will cause the bulge component to be overesti- 
mated. The effects of clumps or asymmetric features on 
the extracted smooth 2D profile parameters were studied 
by adding an asymmetric light component, in the form of 
one or multiple "blobs or HII regions," i. e., unresolved 
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sources convolved with the PSF, to simulated smooth 2D 
profile images. 

The input parameters for generating the asymmetric 
features are the number of HII regions tihii , the total flux 
in the HII regions as a fraction fan of the total galaxy 
flux, and the HII regions' maximum galactocentric dis- 
tance Thii in units of galaxy half-light radius. The posi- 
tions of the HII regions were randomly distributed within 
a circular aperture defined by run- No overall clliptic- 
ity or radial exponential weight was given to the spatial 
distributions of the HII regions. A radially weighted el- 
liptical HII region distribution could originate with HII 
regions linked with an inclined galaxy disk, but the "HII 
regions" here are meant to represent all unresolved asym- 
metric structures, and asymmetric structures may or may 
not necessarily be associated with the disk components of 
galaxies. 

The total asymmetric flux was distributed among HII 
regions using a simple recipe (below) with no attempt to 
include, say, a realistic HII luminosity function. The sim- 
ple recipe produced asymmetric structures that visually 
looked reasonable. According to this adopted recipe, the 
flux Fmi,i allocated to the i th HII region was generated at 
random between and 



and size offsets in each magnitude-size range are in agree- 
ment with the systematic errors shown in Figures 4 and 5. 
The bulge fraction is also fairly robust against asymme- 
tries. The median bulge fraction error is only about 0.1- 
0.2 for galaxies with fun = 0.20 — 0.25. However, in- 
creasing asymmetric flux leads to increasing scatter in the 
bulge fractions, and this scatter is skewed towards lower 
bulge fractions i. e., bulge fraction is always underesti- 
mated when asymmetries matter in a galaxy. The asym- 
metry parameter Ra + Rt (bottom half of Figure 14) re- 
covers most of the asymmetric flux in big, bright galax- 
ies. Equally bright but smaller galaxies have measured 
Ra + Rt slightly lower than big galaxies possibly due to 
the fact that the centroid of the bulge+disk models was 
allowed to vary by ± 1 pixel (±0"15) in the fits, and a 
shift in centroid would always be used by the fitting algo- 
rithm to reduce the overall amount of asymmetry "seen" 
by the smooth model. The scatter in the recovered val- 
ues of Ra + Rt is higher for higher asymmetric fluxes due 
to model centroiding errors introduced by the asymme- 
tries themselves. Ra + Rt increasingly underestimates the 
asymmetric flux at fainter and fainter magnitudes as indi- 
vidual asymmetry sources become too faint to be picked 
out of the noise. 
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where 1 < i < Uhii, ^totai.gaiaxy is the total model galaxy 
flux, -Fhij^ is the maximum flux available to the i th HII 
region, Fuuj is the flux that was actually allocated to the 
jth jjjj region, and F H ii j= o = 0. So, the bracket in Equa- 
tion 21 contains the total unallocated HII flux that remains 
after i — 1 regions have been created. For i = nun, ^hii,» 
is automatically set to the left-over HII flux. 

Asymmetric features superposed on the smooth profile 
were generated randomly for Hhii = 5 and trii = 1-5^/. 
Fifteen identical models were created for each simulated 
galaxy. The first five models had no HII regions, and they 
were used to establish a comparison baseline. The remain- 
ing ten models were divided into five discrete flux levels 
Uhii = 0.05, 0.10, 0.15, 0.20, 0.25) with two models at 
each level in order to sample the same range of residual 
fluxes as seen in the real GSS galaxies. This set of sim- 
ulations contains 170 different galaxy models for a total 
of 2550 simulations. The galaxy models were created with 
the following structural parameters: m^8i4iy (AS) = 24.0, 
B/T = 0.3, r e = 0"12, e = 0.2, r d = 0:'32, i = 20°, 
4>b = 4>d = 60°, and n = 4.0. These asymmetric image 
galaxy simulations were analyzed exactly the same way as 
the real data, and the biases in the parameter values re- 
covered by GIM2D were then examined at each fun level. 

Figures 13 and 14 show the median systematic error on 
the recovered B/T, ru, and Jgi4 parameters as function 
of fmi in different six magnitude-size ranges. Both the 
recovered total magnitude Igu and half-light radius r^i 
show no trends with increasing fan, and the magnitude 



8. SURVEY SELECTION FUNCTIONS 

Generalizing the formalism developed in Simard et al. 
(1999), the observed distribution of galaxies in structural 
parameter space as a function of redshift ^o(W, z) is the 
result of any inherent changes in the resident 11 galaxy dis- 
tribution ^[/(W, z ) in that space and of observational se- 
lection effects. Selection effects are likely to be significant 
given the wide range of structural parameters observed 
locally (Bender, Burstein & Faber 1992; Burstein et al. 
1997). It is therefore important to carefully characterize 
selection effects to disentangle them from real changes in 
<J>i/(W, z). The path from ^u(W, z) to *o(W, z) is given 
by: 

*o(W, z) - Sps(W, z)S UP (W, z)v%(W, z), (22) 

where W is the full set of intrinsic structural parame- 
ters (note the use of lower and upper cases to distinguish 
between apparent and intrinsic structural parameter sets 
here). The subscript UP stands for "Universe to Photo- 
metric sample," and the subscript PS stands for "Pho- 
tometric sample to Spectroscopic sample." The resident 
galaxy distribution ^fu(W, z) is not known a priori. Once 
the two selection functions in Equation 22 have been char- 
acterized, their product (denoted Sus hereafter) shows 
the volume of the structural parameter space where real 
galaxies would have been observed if they existed in that 
region at high redshift. The spectroscopic selection func- 
tion Sps(W, z) is derived in Phillips et al. (2002), so the 
remainder of this section will focus on Sup(W, z). 

The selection function Sup(W, z) contains the informa- 
tion needed to go from any sample of galaxies on the sky 
to the photometric catalog produced with SExtractor and 
reflects the adopted SExtractor detection parameters (de- 
tection threshold in sigmas, minimum detection area, etc.). 



11 It is very important to note the use of the term "resident" here and throughout the rest of the paper to refer to the intrinsic galaxy population 
at a given redshift z. In the absence of real evolution in the galaxy population with redshift, all resident populations would be the same as the 
local population of galaxies. 
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The detection thresholding method used by SExtractor 
depends critically on galaxy apparent surface brightness. 
The probability that a given object will be detected de- 
pends on total flux F, bulge fraction B/T, photobulge ef- 
fective radius r e , photobulge ellipticity e, photodisk scale 
length rd, and photodisk inclination i. For example, ob- 
jects with larger B/T will be easier to detect because they 
are more concentrated, and large objects will be harder 
to detect than smaller ones at a fixed total flux. The se- 
lection function does not depend on the photobulge and 
photodisk position angles. However, note that the selec- 
tion function will also depend on disk internal extinction 
(if any), but this dependence is neglected here since GSS 
galaxies were analyzed with optically thin disks (C a bs = 0, 
recall Equation 7 in Section 5.2). In practice, Sup(W , z) 
is derived from the selection function Sup(w) determined 
as a function of the observed structural parameters. The 
transformation Sup(w) — * Sjjp(w, z) can be made in each 
rcdshift bin using fc-corrections calculated with the median 
observed galaxy Ve;o6 — hit color of the B/T < 0.2 galax- 
ies at that redshift and the cosmological scale relations for 
the assumed cosmology. 

Sup(w) was constructed by generating 60,000 galaxy 
models with structural parameter values uniformly cover- 
ing the ranges: 20.0 < J 8 i4 < 25.0, 0.0 < B/T < 1.0, 
0"0 < r d < 10"0, < sin i < 0.9962. It is better 
here to uniformly generate disk inclinations in sin i rather 
than in i since randomly oriented, optically thin disks in 
space are expected to have a uniform sin i distribution. 
Each model galaxy was added, one at a time, to an empty 
20" x 20" section of a F8UW HST/GSS image (same im- 
age section as used in Section 7.1). "Empty" here means 
that no objects were detected by SExtractor in that sky 
section with the same detection parameters used to con- 
struct the object catalog. Using an empty section of the 
GSS ensured that Sjjp(\v) was constructed with the real 
background noise that was seen by the detection algo- 
rithm. The background noise included read-out, sky and 
the brightness fluctuations of very faint galaxies below the 
detection threshold. This last contribution to the back- 
ground noise is particularly hard to model theoretically, 
and the current approach bypassed this problem. SEx- 
tractor was run on each simulation with the same param- 
eters that were used to build the SExtractor catalog. The 
function Sup(w) was taken to be the fraction of galaxies 
successfully detected and measured by SExtractor at each 
location w in structural parameter space. 

Figure 15 shows one-dimensional projections of S[/p(w) 
onto each of the six structural parameters (F, B/T, r e , 
e, rd, and i). Different symbols show groups of galaxies 
with different bulge fractions. Clearly, some parameters 
are more important for the selection function than oth- 
ers, and bulge parameters are obviously more important 
for bulge-dominated galaxies and the same holds true for 
disk parameters. The selection function depends strongly 
on total apparent magnitude independent of bulge frac- 
tion, and bulge-dominated galaxies are more likely to be 
detected at all magnitudes than disk galaxies. The selec- 
tion function for bulge-dominated galaxies decreases with 
increasing bulge effective radius out to r e = 2" and re- 
mains relatively flat beyond that radius. There is a very 
weak dependence of Sup(e) on bulge ellipticity. The se- 



lection function for B/T > 0.2 galaxies is nearly indepen- 
dent of disk scale length whereas the selection function for 
B/T < 0.2 galaxies decreases out to r& = 2.2" and re- 
mains flat after that. There is virtually no dependence of 
Sup on the disk inclination angle, and this is somewhat 
surprising given that more inclined, optically thin disks 
should be easier to detect. This apparent puzzle was re- 
solved by generating a set of 60,000 pure disk galaxy mod- 
els and re-computing Sup with this new set. There was a 
clear dependence of Sup on disk inclination for this pure- 
disk galaxy set. The detectability of disks went from 0.55 
at i ~ 0° to 0.75 at i ~ 80°. Given that the disk sample in 
the selection function shown in Figure 15 includes galaxies 
with bulge fractions between 0.0 and 0.2, it appears that 
even a small (luminosity-wise) de Vaucouleurs bulge can 
boost the detectability of a galaxy enough to mask out the 
disk inclination dependence of the selection function. 

9. COMPARISON WITH THE MEDIUM DEEP SURVEY 

The Medium Deep Survey (MDS, Ratnatunga et al. 
1999 and references therein) is the largest database of HST 
galaxy structural parameters in existence with 200,000 ob- 
jects as of October 1998. The images analyzed by the 
MDS team consist of MDS WFPC2 pure parallel observa- 
tions as well as of HST archival observations of randomly 
selected WFPC2 fields such as the Groth Strip and the 
Hubble Deep Field among others. The MDS team fitted 
the profiles of Groth Strip galaxies separately in Veoe and 
/8i4 using a completely different analysis pipeline, a com- 
pletely different likelihood maximization algorithm and a 
different bulgc+disk model. The parameters of the MDS 
bulgc+disk galaxy model are sky background, x-y cen- 
troid, orientation (bulge and disk are assumed to have the 
same position angle), bulge and disk axis ratios, bulge frac- 
tion B/T and the ratio of the bulge/disk half-light radii. 
The large sizes (thousands of objects each) of the MDS and 
DEEP/GIM2D Groth Strip structural catalogs make them 
ideally suited to run a check of one against the other. This 
is the first time that MDS results are compared against an 
independent work on such a scale. 

The MDS Maximum-Likelihood Estimate (MLE) struc- 
tural catalogs of the Groth Strip galaxies were extracted 
directly from the on-line MDS CD-ROMs, and the MDS 
catalogs were matched to the DEEP/GIM2D separate 
structural fit catalog using a matching radius of 0"8 and 
a maximum 7si4 magnitude difference of 1. The match 
yielded 7138 positive cross-identifications. The results of 
the match are shown in Figure 16. The top left-hand panel 
shows the GIM2D galaxy model Isu total magnitudes 
against the MDS galaxy model Isi4 total magnitudes. The 
long-dashed line is a one-to-one line, the filled circles are 
galaxies with \I S u,gim2D - hu,MDs\ < 0.2 mag, and the 
open circles are galaxies with magnitude differences larger 
than 0.2 mag. The envelope of the data point distribution 
is clearly asymmetric with respect with the one-to-one line 
in its upper section i. e., GIM2D magnitudes for some ob- 
jects are too faint with respect to MDS magnitudes. The 
asymmetry in the upper envelope is due to the fact that it 
is made up of two distributions. One distribution comes 
from real photometric errors, and the distribution of its 
points as a function of distance from the one-to-one line is 
symmetric with respect to the lower envelope. The second 
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distribution inside the upper envelope contains the most 
discrepant objects, and it is due to the fact that objects 
are more finely split in the DEEP/GIM2D structural cat- 
alog than in the MDS catalog. Thus some single objects 
in the MDS catalog are two or more distinct objects in the 
DEEP / GIM2D catalog, and their DEEP/GIM2D magni- 
tudes are therefore fainter. 

The top right-hand panel of Figure 16 shows the GIM2D 
galaxy log half-light radii in arcseconds against the MDS 
radii for galaxies with I^u < 22. The long-dashed line is 
again the one-to-one line. The half-light radii are in excel- 
lent agreement with no systematic differences. The mean 
log radius difference (GIM2D— MDS) and rms scatter are 
(0.002, 0.049) for objects with AJ 8 i 4 < 0.2 mag (filled 
circles) and (-0.071, 0.168) for objects with A7 8 i 4 > 0.2 
mag (open circles). The rms scatter for the open circles 
is higher as one would expect from using different object 
splitting. 

The GIM2D /si4 bulge fractions are compared against 
MDS Igi4 bulge fractions for galaxies with /gi4 < 22 in 
the last panel of Figure 16. The two vertical distribu- 
tions at (B/T)mds — 0.0 (pure disk systems) and at 
(B /T)mds = 1-0 come from the MDS pipeline where ob- 
jects below a certain S/N threshold and/or size are only fit- 
ted by a pure bulge or a pure disk model and not by the full 
bulge+disk model. Setting these points aside for now, the 
mean B/T difference (GIM2D-MDS) and rms scatter for 
intermediate B/T galaxies are (—0.021,0.105) for the filled 
circles and (—0.046, 0.235) for the open circles. The slight 
systematic offset and the increased rms scatter for the open 
circles is not surprising since objects that were split dif- 
ferently are likely to be classified differently. The discrep- 
ancies between the MDS and GIM2D bulge fractions at 
{B/T) mds — and {B/T)mds — 1 ar e best studied by 
looking at the two extremes where (B/T)gim2D > 0.5 at 
(B/T) MDS = and {B/T) GIM2D < 0.5 at {B/T) M ds = 
1. In the first case, the objects are unresolved in both MDS 
and GIM2D catalogs. In the second case, there are differ- 
ent reasons behind the bulge fraction discrepancies. Nine- 
teen objects have (B/T)gim2D < 0.5 at (B/T)mds = 1, 
and each was visually inspected. Out of these 19 objects, 
6 objects are irregular/peculiar galaxies or close pairs not 
separated in the MDS catalog, 2 are stars, 5 are single ob- 
jects with (B/T) MDS jsi4 = 1 and (B /T) M ds,V606 < 0.5, 
3 are equally well fit by the MDS and GIM2D models, 2 
are single objects with knots in their MDS Isu residual 
images, and 1 is a bad object match. Even though these 
discrepant objects are very interesting, they represent a 
very small fraction of both catalogs, and the MDS and 
GIM2D results are generally in good agreement. 

10. SUMMARY 

The structural parameters of galaxies in the Groth Sur- 
vey Strip were measured from archival HST images as part 
of a combined HST and Kcck/LRIS survey of the Strip by 
the DEEP team. GSS galaxy surface brightness distri- 
butions were fitted with a 2D, PSF-convolved bulge+disk 
model using an implementation of the Metropolis algo- 
rithm to optimize model parameters. A total of 7450 
galaxies were fitted separately in 7 8 i4 and Vqo6 ■ 648 galax- 
ies with secure Keck/LRIS redshifts were also fitted simul- 
taneously in both bandpasses. The structural catalogs in- 



clude image asymmetry parameters and rest-frame magni- 
tudes and colors for bulge and disk components computed 
using two different sets of fc-corrections. 

This paper also provides full structural catalogs of two 
extensive sets of close to 6000 simulations to allow de- 
tailed characterizations of the systematic and random er- 
rors in separate and simultaneous structural fits at any 
location in structural parameter space. Error maps for 
galaxy magnitudes, half-light radii and bulge fractions are 
presented as examples. The simultaneous structural fit 
simulations include varying bulge and disk colors and show 
that bulge and disk colors can be reliably measured down 
to Isu(bulge) = 23.5 and Isu(disk) — 23.5. Similar for- 
mats are used for the real and the simulation catalogs so 
that interested users can study biases associated with dif- 
ferent selection criteria to see which criteria will make the 
best use of the real data for the specific science goals being 
pursued. 

The effects of unresolved, "HII-region"-like asymmet- 
ric structures on fitting parameters were studied with a 
third set of 2550 simulations. Recovered total magnitudes, 
half-light radii and bulge fraction were, on average, robust 
against the presence of strong asymmetries with little or 
no systematic bias. However, the scatter in the recovered 
parameters increased with increasing asymmetric flux, and 
this increased scatter was skewed towards lower bulge frac- 
tions. Bulge fractions are always underestimated when the 
presence of strong asymmetries matters. The asymmetry 
parameter Ra + Rt was found to be a good estimate of 
the total asymmetric flux present in large, bright galaxies. 

The photometric selection function of the survey was 
mapped over a wide range of magnitudes, bulge fractions, 
and bulge/disk sizes to delineate the volume of structural 
space favored by the source detection algorithm. The 
key quantity is surface brightness. At a given magnitude, 
larger galaxies are harder to detect, and disks are harder 
to detect than bulges since disk profiles are less compact. 
There is little or no dependence of the selection function on 
bulge cllipticity and disk inclination. The selection func- 
tion, coupled with biases and errors from the simulations, 
gives a complete reproduction of the observational limits. 

The structural parameters presented here were com- 
pared with the results from the Medium Deep Survey 
database. This is the first large scale comparison of 
MDS results against an independent source. The MDS 
and DEEP/GIM2D catalogs yields 7138 positive cross- 
identifications, and measurements of total magnitude, 
half-light radius and bulge fraction are all in good agree- 
ment with no systematic differences. Bulge fraction clas- 
sifications are in disagreement for a small fraction of the 
galaxies. 

The combination of the three essential ingredients of 
quantitative galaxy morphology (very large sets of struc- 
tural measurements, detailed characterization of biases 
and errors, and mapping of the multidimensional photo- 
metric selection function) found in this paper is ideal to 
pursue direct comparisons between observations and the 
latest models of galaxy formation and evolution. 

The development of GIM2D greatly benefited from the 
experiences of a core user group which includes Kim-Vy 
Tran, Brad Holdcn, Francine Marleau, Kathcrine Wu, 
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Fig. 1.— HST/WFPC2/F8UW image of the Groth Survey Strip Field 8/Wide Field Camera Chip 3. Image scale is O'.'l pixel" 1 , and 
the total exposure time is 4400 seconds. 

Fig. 2. — SExtractor segmentation image of CSS image in Figure 1. Every pixel in this image was assigned the non-zero SExtractor flag 
value of its parent object or a flag value of zero if it was a background pixel. GIM2D uses this segmentation image to independently fit objects 
which are very close to each other. 

Fig. 3. — Mosaics of CSS thumbnail images drawn at random from the DEEP/GSS sample. The extracted region for each galaxy is twenty 
times the area of its 1.5— <r bkg isophote. First mosaic - Science thumbnail images: The numbers shown around each science thumbnail image 
are the DEEP/GSS object ID (top left), DEEP Keck redshift (top right), observed galaxy total magnitude /gi4 (bottom left) and observed 
galaxy Veoe — -?814 color (bottom right). Second mosaic - Output PSF-convolved GIM2D model thumbnail images: The numbers shown around 
each model thumbnail image are the DEEP/GSS object ID (top left), bulge fraction B/T (top right), rest-frame B-band galaxy total absolute 
magnitude (bottom left), and galaxy log half-light radius in arcseconds (bottom right). Third mosaic - Residual thumbnail images (science 
— model): Numbers shown around each residual thumbnail image are the DEEP/GSS ID (top left), RT3 + RA3 asymmetry parameter (top 
right), concentration index C (bottom left) and asymmetry index A (bottom right). The same greyscale was used for the science, model and 
residual thumbnail images of a given galaxy, but a different greyscale was used for each galaxy. The bottom greyscale value for each galaxy 
was set to be 5-ai,kg below the science thumbnail background level, and the top greyscale value was set to be 10-<T(,fc 9 above the background. 
Notice the close galaxy pair 192.2330 which was successfully deblendcd by SEXtractor and independently fitted by GIM2D. Also notice the 
wealth of interesting structures in all the residual images. 

Fig. 4. — Two-dimensional maps of systematic and random galaxy magnitude errors from 5995 DEEP/GSS separate structural fit sim- 
ulations. Top left-hand panel: Systematic error on recovered galaxy total magnitude Isi4,rec as a function of input galaxy log half-light 
radius r^j iinpllt in arcseconds and input galaxy total magnitude ^814, input- The top number in each cell is the mean magnitude error 
(^814, rec — ^814, input), and the bottom number is the number of simulations created in that cell. Top right-hand panel: l-cr random error 

on 7814, rec {&{Is\4,rec — ^814, input)) 3- function of log T hi, input 

and ^814, input- Bottom left-hand panel: Systematic error on 7814, rec as a 
function of recovered galaxy log half-light radius rj,j ]7 . ec in arcseconds and Isi4,rec- Bottom right-hand panel: Actual DEEP/GSS /8l4,o6s-l°g 
r hl,obs observations. Since recovered quantities are the simulation equivalent of the observed ones, this panel should be directly compared to 
the bottom left-hand panel to see how important magnitude errors are in different regions of the observational magnitude-size plane. Although 
such 2D maps can be very useful tools, one should keep in mind that they are only two-dimensional projections of a complex multi-dimensional 
error function. 

Fig. 5. — Same as Figure 4 except that log half-light radius errors are shown here. 
Fig. 6. — Same as Figure 4 except that bulge fraction errors are shown here. 

Fig. 7. — Two-dimensional maps of systematic and random galaxy bulge fraction errors from DEEP/GSS separate structural fit simulations 
with input galaxy half-light radius r^i^nput > O'.'l. Top left-hand panel: Systematic error on recovered galaxy bulge fraction (B/T) rec as a 
function of input galaxy bulge fraction (B /T)i nput and input galaxy total magnitude Igi4,input- The top number in each cell is the mean 
bulge fraction error ((B /T) rec — (B/T)i npu t), and the bottom number is the number of simulations created in that cell. Top right-hand panel: 
l-cr random error on (B/T) rec (cr((B/T) rec — (B/T)i„ pu t)) as a function of (B /T)i npu t and ^814, input- Bottom left-hand panel: Systematic 
error on (B/T) rec as a function of (B/T) r ec and Isi4,rec- Bottom right-hand panel: Actual DEEP/GSS (B/T) j, s -/8l4,o6s observations. 
Since recovered quantities are the simulation equivalent of the observed ones, this panel should be directly compared to the bottom left-hand 
panel to see how important bulge fraction errors are in different regions of the observational magnitude-bulge fraction plane. 

Fig. 8. — Two-dimensional maps of systematic and random galaxy magnitude errors from 5195 DEEP/GSS simultaneous structural fit 
simulations. Top left-hand panel: Systematic error on recovered galaxy total magnitude 7814, rec as a function of input galaxy log half- 
light radius r^j input m arcseconds and input galaxy total magnitude Igi4, input- The top number in each cell is the mean magnitude error 
(^814, rec — Isi4,input), and the bottom number is the number of simulations created in that cell. Top right-hand panel: l-cr random error 
on ^814, rec (o~(Isi4,rec — I&i4,input)) as a function of log r^i input and /814, input- Bottom left-hand panel: Systematic error on ^814, rec as a 
function of recovered galaxy log half-light radius r^i^rcc in arcseconds and i8l4,rec- Bottom right-hand panel: Actual DEEP/GSS /8i4,o6s-log 
r hl obs observations. The number of data points is smaller than in Figures 4, 5, and 6 since simultaneous bandpass fits were performed only 
on GSS galaxies with secure DEEP Keck rcdshifts. Since recovered quantities are the simulation equivalent of the observed ones, this panel 
should be directly compared to the bottom left-hand panel to see how important magnitude errors are in different regions of the observational 
magnitude-size plane. Although such 2D maps can be very useful tools, one should keep in mind that they are only two-dimensional projections 
of a complex multi-dimensional error function. 

Fig. 9. — Same as Figure 8 except that log half-light radius errors are shown here. 

Fig. 10. — Same as Figure 8 except that bulge fraction errors are shown here. 

Fig. 11. — Structural component colors recovered by GIM2D from simultaneous structural fit simulations. Top three panels: Recovered galaxy 
^606 ~ Isi4 color versus input galaxy Vqo6 ~ ^814 color for input galaxy total magnitudes in the ranges 20.0 < isi4 < 21.0, 21.0 < /si4 < 22.0, 
and 22.0 < / 8 14 < 23.5. The mean and rms color differences (recovered-input) are (0.002, 0.020), (0.010, 0.034), and (0.063, 0.087) for the 
three ranges respectively. Vgo6 lim = 26.0. Middle three panels: Recovered bulge Vqo6 — 1&14 color versus input bulge Vgog — Igi4 color for input 
bulge magnitudes in the ranges 20.0 < /8i4(bulge) < 21.0, 21.0 < /8l4(bulge) < 22.0, and 22.0 < 7si4 (bulge) < 23.5. The mean and rms bulge 
color differences (recovered— input) are (-0.001, 0.255), (0.022, 0.120), and (0.034, 0.186) for the three ranges respectively. Vqqq Um (bulge) = 
26.0, r^i gi4 > 0"15, and r^i qqq > 0"15. Bottom three panels: Recovered disk Vaoe — Isi4 color versus input disk Vqog ~ I&14 color for input 
disk magnitudes in the ranges 20.0 < isi4(disk) < 21.0, 21.0 < i8l4(disk) < 22.0, and 22.0 < /8l4(disk) < 23.5. The mean and rms disk color 
differences (recovered— input) are (0.008, 0.047), (0.006, 0.211), and (0.030, 0.304) for the three ranges respectively. V606, lim (disk) = 26.0, 
»*M,814 > 0"15, and r w ,606 > 0"15. 

Fig. 12. — Two-dimensional maps of systematic galaxy bulge fraction errors from DEEP/GSS separate and simultaneous structural fit 
simulations as a function of input galaxy bulge fraction and input bulge effective radius/disk scale length ratio. Left-hand panel: Systematic 
error on recovered galaxy bulge fraction (B/T) rec from separate structural fits (SPF) as a function of input galaxy bulge fraction (B/T)i npu t 
and input log bulge effective radius/disk scale length ratio re/r^. The top number in each cell is the mean bulge fraction error ((B /T) rec — 
(B /T)i npu t), the middle number in each cell is the minimum error, and the bottom number in each cell is the maximum error. Right-hand 
panel: Same as left-hand panel but for simultaneous structural fits. These maps show that there do not seem to be any regions where GIM2D 
can systematically mistake disks for bulges and vice versa. 

Fig. 13. — Median systematic errors on recovered bulge fraction B/T and galaxy total apparent magnitude 7si4 as a function of asymmetric 
flux fraction full- Top six panels: Median B/T systematic error in six different magnitude-size ranges. Bottom six panels: Median Isi4 
systematic error in six different magnitude-size ranges. Vertical error bars are 68% lower and upper bounds. 
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Fig. 14. — Top six panels: Median systematic error on recovered galaxy log half-light radius r^; as a function of asymmetric flux fraction 
fmi in six different magnitude-size ranges. Vertical error bars are 68% lower and upper bounds. Bottom six panels: Median recovered 
Ra + Rt asymmetry index values as a function of asymmetric flux fraction full in six different magnitude-size ranges. Vertical error bars 
arc 68% lower and upper bounds. 

Fig. 15. — Six one-dimensional projections of the photometric selection function Sjjp(v7): galaxy total apparent magnitude 7gi4 (top left), 
galaxy bulge fraction (B/T) (top middle), bulge effective radius r e in arcseconds (top right), bulge ellipticity e (bottom left), disk scale length 
in arcseconds (bottom middle), and disk inclination angle i (bottom right). Filled circles are simulated galaxies with 0.0 < (B/T) < 0.2, 
pluses are galaxies with 0.2 < (B/T) < 0.8, and triangles arc galaxies with 0.8 < (B/T) < 1.0. 

Fig. 16. — Comparison between MDS and GIM2D parameters for Groth Strip galaxies. Top left: GIM2D versus MDS total galaxy model 
7gl4 magnitudes. Top right: GIM2D versus MDS galaxy semi-major axis log half-light radiii in arcseconds for galaxies with igl4 < 22. Lower 
left: GIM2D versus MDS galaxy bulge fractions for galaxies with Isu < 22. Filled circles are galaxies with 1^814, GIM2D — Isu,mds\ < 0.2, 
and open circles are the remainder of the sample. The long dashed lines in all three panels are one-to-one lines. 
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